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Cet ouvrage propose une collection d’algorithmes numériques, destinés à 
être employés sur micro-ordinateur dans le but de résoudre nombre de 
problèmes scientifiques, relatifs aux mathématiques, à la physique et à 
l'électronique, et s’adresse à tous les possesseurs de micro-ordinateurs, 
quel que soit leur type. Tous ces programmes ont d’ailleurs été mis au 
point sur une version 16 K. Ils sont écrits en BASIC standard et ne font 
pas appel aux propriétés particulières d’une machine. 

Ce livre s'adresse à tous les curieux, bloqués devant un problème qu'ils 
ne savent pas résoudre. C’est alors que l’on s’aperçoit que l'explosion de 
l'informatique, et les nouveaux modes de pensée qu’elle devrait induire, 
n'ont pas toujours été accompagnés par la réforme nécessaire de 
l’enseignement théorique dans nos lycées et facultés. 

Prenons l'exemple (ô combien classique!) du pendule pesant, au 
programme des classes de terminale. Les connaissances des élèves leur 
permettent aisément d’écrire l’équation différentielle du mouvement; il 
s’agit en effet d’un problème simple de mécanique. 

Un double écueil surgit alors : d’une part, cette équation n'est pas 
linéaire — on décide alors, assez arbitrairement, de n'étudier que les 
«petits» mouvements —; d'autre part, les connaissances mathématiques 
des élèves permettant de résoudre cette nouvelle équation linéarisée ne 
sont pas du niveau de cette classe; les mathématiciens fournissent alors la 
solution (signalons à leur décharge qu’ils peuvent mener une discussion 
qualitative de l'équation exacte et décrire les conditions initiales 
conduisant à des solutions périodiques). 

Or, grâce à l'informatique, la seule connaissance de la première 
équation, établie par les élèves et comprise par eux, permet d'aboutir, 
sans aucune approximation, à la solution (numérique), seul objet digne 
d'intérêt pour le physicien. 

Qu'importe la méthode! Et n'est-il pas merveilleux de trouver en 
quelques secondes la solution de problèmes que des dizaines de 
générations ont étudiés péniblement? 

Nous nous sommes particulièrement attachés à proposer des program- 
mes courts et modulaires, sans pour autant sacrifier la précision ou la 


AE 


rigueur des algorithmes; nous expliquons toujours largement la méthode 
employée. En outre, quand cela était possible, nous avons éliminé les 
recherches de convergence (la solution existe-t-elle?), puisque tout 
problème physique a une solution. Tous les programmes ont naturelle- 
ment été testés, notamment sur le choix d'exemples proposés, et sont 
utilisables directement pour d’autres applications. Aucune digression 
théorique n'est faite, surtout si elle ne débouche que sur une méthode 
inexploitable. Un dernier problème, pratique, subsiste : toutes les 
méthodes proposées mettent d'autant plus de temps à s’effectuer qu’une 
grande précision leur est demandée. Il est alors conseillé de travailler sur 
un micro-ordinateur disposant d’un BASIC compilé, ou, si ce n’est pas 
possible, d’un Pascal. Nous donnons d’ailleurs les programmes les plus 
«lents» dans ces deux langages. 

Enfin, tous ces programmes sont interactifs et se complètent. Nous 
donnons à titre d'exemple la combinaison des Chapitres 5 et 6, qui 
permet d'obtenir d’une façon très élégante la solution d’un problème 
particulier. 

Nous remercions nos élèves respectifs de nous avoir incités à écrire ce 
document et de nous avoir guidés dans sa réalisation. 


PIERRE-MARC BEAUFILS, 
WOLFRAM LUTHER. 


Les instructions relatives au mode texte sont standard. Pour le mode 
haute résolution, les caractéristiques sont les suivantes (elles demandent 
donc éventuellement une adaptation) : 


Axe OX horizontal gradué de 0 à 239. 

Axe OY vertical gradué de 0 à 199. 

Origine des axes située en haut à gauche. 

CURSET X, Y, 1 : allume un pixel de coordonnées absolues X et Y. 


CURMOV X, Y, 1 : déplace le point courant des quantités X et Y 
relativement à l’ancienne position. 


e DRAW X, Y, 1 : trace une droite de coordonnées X et Y à partir du 
dernier point. 


e CHAR A, 0, 1 : écrit le caractère de code ASCII A à la position du 
point courant. 


e CIRCLE R, 1 : trace un cercle de rayon R ayant pour centre le dernier 
pixel allumé. 


| 
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1. NOMBRES COMPLEXES. MODULE ET PHASE 





En électronique, il est souvent fait usage de la notation complexe pour 
représenter une grandeur relative à un fonctionnement d’un système en 
régime permanent sinusoïdal. En particulier, la transmittance d’un 
quadripôle, rapport entre tension de sortie et tension d’entrée, est une 
grandeur complexe, fonction d’un paramètre, la pulsation du signal 
sinusoidal incident. De façon générale, et dans les cas les plus simples, 
cette transmittance peut être de trois types fondamentaux : 


1. passe bas du second ordre (1), 
2. passe haut du second ordre (2), 
3. passe bande du second ordre (3). 


Leur expression mathématique est la suivante : 





T, = à 
. © 
1+2 jm 24 (52) 
0 (O7) 
; (0) 
2 jm — 
— (OR 1 
RE 
.. & . © ] DO 
1+2im + (52) 14 (2-2) 
0 6 2m\a © 


Les électroniciens étudient ces fonctions dans le plan de Bode, plan 
dans lequel l’axe des abscisses représente le logarithme de la pulsation et 
celui des ordonnées 20 log|T| (en dB donc). 

Tracer ces fonctions ne présente pas de difficultés particulières. 
Cependant, l’étude de la phase pose un problème dans la mesure où le 
dénominateur commun à ces transmittances a un argument variant de 0 à 
— 7 lorsque &« varie de O0 à l'infini. Pour la pulsation w,, on constate en 
effet que la tangente de cet argument devient infinie. Il est alors 
préférable de se servir de la seconde forme de l’expression de T, qui, elle, 
ne possède pas de singularité gênante, et d'écrire : 


— T — 
Arg(T;)= nc Arg(T:) 


== T = 
Arg(T;) = 3 + Arg(T;) 


Nous proposons de tracer les lieux de Bode de ces trois types de filtre. 
Il est possible d’en déduire la réalisation de courbes universelles. 

La construction de l’échelle logarithmique est effectuée de la manière 
suivante : 

Supposons que 50 points représentent une décade. Nous définissons un 


1 
incrément D= 107 (=) = 10 0,02. Nous passons donc d’un point au 


suivant en multipliant la fréquence par cet incrément. Il paraît ainsi 
évident qu’au bout de 50 points : f'=f.(D)Î 50—=10f. On peut donc 


1 
dire que la distance séparant deux points représente 30 de décade. 


Nous proposons dans le programme deux graduations possibles pour 
cet axe horizontal, soit 100 points = 1 décade (on observe ainsi l’évo- 
lution des courbes sur deux décades), soit 50 points/décade (observation 
de quatre décades). 
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14 CLS 
£4 PRIHT'TYPE OÙ FILTRE :" 


34 PRINT 
49 PRINT" #Passe bas... ,.,,612" 
59 PRINT" #Passe haut. ..,.,625" 
64 PRINT" #Passe bande... . 6x5" 
65 GET KR 
78 PRINT 


84 INPUT'"'Facteur d'amortissement: ":t 
94 IF MSITHENTR=ZAKLOGÉ 2#M/SORC 1-28 5 à 


ELSE TR=A 

92 PRINT:PRIHT'"'Voulez vous etudier la f 
onction sur d decades( 15 ou 2 decadesé(2 
7": GETH 


199 IF H=1 THEN D=19"9,42 ELSE D=19"a.a 
1 


PE IF N=1 THEH W=1/-<19440% ELSE W=1/c1 
G#D 

124 ON K GOSUR 344,3274,34g 

139 DEF FHPCXUS-ATME Ce 8 176 SEM 5 +P 
149 FRINT:PRINT'"'Voulez vous le module(h 
3 ou 14 Phase PiT" GET AS:HIRES 

145 GOSUB 1964 

159 IF A$S="M'TMHEH 169 ELSE IF A$="P" TH 
EN 299 

153 GOTO 148 

169 FOR *“=9 TO 193 

179 CURSET20+%, 60-2A#LO0GC FMMCU 5 5, 1 

18 W=W#D 

139 NEXT 

135 STOP 

284 FOR X=A TO 199 

219 CLIRSET29+%, 60-< 30/€ PI /2 5 3#FNFCU 5, 1 


228 W=H#D 

234 NEXT 

244 STOF 

3908 DEFFHME 4 = 4e CC 1-44 10240 ZKMEK 292 3% 
9.53:P=-PI1/2 

319 RETURN 

320 DEF FNMCX m2 CC CL AE DDC 2EMXE 37 
2"0.5):P=P]1/2 
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334 RETURN 

344 DEF FMC KO Z2#MAK 6 CC LH DDC DAMES 
»22"9.5):Peg 

350 RETURN 

1909 REM ##%# AXES #X* 

1919 CURSETO, 0,1 :DRAUZ3S, 0, 1:CFAUG, 199. 
1:DRAU-239,0,1:DRAUG,-199,1 

1 29 CURSET20,60,0:DRAU218,09,1:DRAU-3,- 
3:2:CHAR6Z,0,1 

1939 CURSET20,199,0:DRAUG,-199, 1 :OKAU-2 
,=4, A: CHARG4, 0,1 

1949 FOR NH = 1 TO 3 

1959 CURSET 2A+50%N,62,G:DRAUA, -4,1 
1969 NEXT 

197Q FOR N = 9 TO 8 

1089 CURSETZA,20+20#H,06:DRAH 4.4.1 

1994 NEXT 

1199 RETURN 






nent menée nero ne 


20 log |T| 
l 40 dB 
20 dB 


Filtre passe bas. Courbe de module. m = 0,2. Étude sur 4 décades. Tx = 8,13 dB. 


= 45.2 





Filtre passe bas. Courbe de module 2 décades. m—0,2. TR = 8,13 dB. 





Filtre passe bas. m—0,2. Courbe de phase. Étude sur 2 décades. 


2.46 = 










20 log ÎT| 
40 dB 


20 dB 


Filtre passe haut. m = 0,1. Courbe de module. Étude sur 2 décades. Ty = 14 dB. 





Filtre passe haut. Courbe de phase. m=0,1. Étude sur 2 décades. 
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20 log # 20m 
- 40 dB 





Filtre passe bande. m—0,05. Étude sur 2 décades. Courbe de module. 





Filtre passe bande. m = 0,05. Étude sur 2 décades. Courbe de phase. 
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2. ÉQUATION DIFFÉRENTIELLE LINÉAIRE DU PREMIER 
ORDRE 








INITIATION AU PROBLÈME 


Il s’agit de tout système physique décrit par une équation du type : 


Tdy(t) 
dt 





+ y(t)=Bx(t) 


Les solutions d’une telle équation sont parfaitement connues pour des 
signaux x(t) simples : échelon, sinusoïde. Dans les autres cas, il faudra 
une intégration numérique. 





y(t)= yo exp (+5 x (u)exp (©) du 


T 


(Voir le programme Intégration.) Cette équation est celle qui régit le 
fonctionnement des circuits dits du premier ordre en électronique : 


L i(t) 


s(t) e(t) R s(t) 
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_:e(—s (8) di(t) 








i(t) = EMILE 
HEC (9=2@ 
dt 
D'où : D'où : 
ds(t) Le L ds(t) _ 
RC + s(1)= e(4) Sn T0 ED 


De façon générale : 
ds(t) 
T——+s(t)=Be(t) 
dt 


T=constante de temps du système (en secondes). 


B = transmittance statique du système. 


MÉTHODES 


Nous proposons de remplacer les dérivées par les rapports d’accroisse- 
ments finis. Ainsi : 


SR 
T—+s=Be 
A: 


Appelons : 
At=intervalle de calcul; 
5, = valeur de s à l'instant k-At; 
e, = valeur de e à l'instant k-At; 
S,_, = valeur de s à l'instant (k—1)At. 
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Ainsi : 


Sx — Ski 
T——— + 5, =Be,. 
À: k k 


T T 
Sk [1 +2] —=Be, are Sx—1- 


Ainsi, la valeur actuelle de s, peut être prédite, si l’on connaît sa valeur 
antérieure s,_, et la valeur actuelle de l'entrée e,. Pour commencer le 
calcul, il faut naturellement une condition initiale, qui est en général 5, 
(ce qui sera par exemple la tension aux bornes du condensateur à 1 — 0). 
On en déduit facilement l’organigramme. 

Remarquons que e =e(kAt) peut être soit une fonction explicite du 
temps : e. =sin(kAt), soit une fonction quelconque définie par un 
tableau (matrice) : e, =E(k). 









Initialisation : 7, At, 5, 
et B. 






















Introduction de e,. Calcul : 5 = 











Pour les connaisseurs, il va de soi que At doit être petit devant T qui 
caractérise la vitesse de réaction du système. Mais qu’en est-il si e est une 
fonction du temps? Nous reviendrons sur cette question plus loin. 


= 90: 


L'erreur systématique à chaque pas est de l'ordre de (At); l'erreur 
globale après sommation est un infiniment petit de l’ordre de At. 


AMÉLIORATION DE LA MÉTHODE 


La méthode précédente est simple. Essayons de la cerner avec plus de 
rigueur. 
Reprenons l'équation différentielle du système : 


ds(t) 
dt 





T +s(t)=Bet(t) 


Intégrons-la entre (k—1)At et kAt: 





T 
(k—1)At dt 


‘A ds(t) at+ | s(D4r=8 | e(t) dt 


Les deux dernières intégrales peuvent s’évaluer avec la méthode du 
trapèze : 


s(t)ouel(t) 





(k—1)At kAt 


Il est facile de voir que l’aire correspondant à l'intégrale [ s(t) dt peut 
être approchée par celle du rectangle hachuré : 


Sk + Ski 
2 


At 


19912 


D'où : 


Sx + Se = 
RS RE un 


2 2 
et donc : 
At B 
T — — — 
2 2 
DRET Sx—1 env 
T+— T+— 
2 2 


Cette relation, légèrement différente de la précédente, fait intervenir 
ek_1- 

Que pouvons-nous penser de la précision de la méthode et du choix 
optimal de At? 

Pour cela, utilisons cette relation dans un cas où la solution de 
l'équation est parfaitement connue : e(t) sinusoïdal. En notation com- 
plexe, on sait que la transmittance de tels systèmes est du type passe bas : 


Posons alors : 

e, = E exp(jwkAt)=E exp(jwkAt) 

s = So exp(jp)exp(jwkAt)=S exp(jwkAt) 

Sx1 = So exp(jpg)exp(jw(k—1)At)=S exp(jw(k—1)At) 
(La solution, en régime permanent, est en effet sinusoïdale.) 


Le développement de l'expression obtenue précédemment donne 
alors : 


B 
a) 
ie (>) 
LÉ At 


2T 
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Si l’on compare les courbes donnant les modules des fonctions T par 


rapport à &, on obtient : 


20 log TI SE &, 
T 
log 






filtre numérique 


Ces deux courbes sont apparemment différentes, cependant, on 
constate qu’elles coïncident pour les pulsations «petites » : 


oAt 
t&g(x) = x; ns 


Cela nous amène à nous interroger sur le sens physique de la méthode 
numérique étudiée. En effet, du point de vue de l’électronicien, il y a eu 
un échantillonnage à intervalle de temps Af, donc avec la pulsation : 


27 
At 


Or (théorème de Shannon), cela n’a de sens que si la plus haute 
pulsation traitée « est inférieure ou égale à la moitié de la pulsation de 


l’échantillonnage ©, : 


2T TT 


1 
<-— © it L— = 
FDP SAT Ar 


2794. 


et donc la condition cherchée : 


AI 


£a 


En pratique, on prendra un coefficient de sécurité. Les résultats 
obtenus par l'intégration numérique ne seront alors en théorie valables 
que si le développement en série de Fourier de e(t) ne comprend pas 


: : fe , T : 
d’harmoniques de pulsations supérieures à ya Pour que cela ait un sens, 
t 
il faudra donc en général avoir : 


re 


Remarque 


Cette méthode n’est plus approchée, mais exacte, si l’on veut simuler 
le fonctionnement d’un filtre numérique, construit suivant la disposition : 





€k 


Sk 







retard 





ligne à 
retard 


At 





Il pourrait sembler que la méthode est d’autant plus précise que le pas 
de calcul est petit. Il n’en est rien. En effet, il n’est pas impossible, si ce 
pas est effectivement petit, qu’un calcul intermédiaire conduise à une 
variation d’une quantité plus petite que la précision du micro-ordinateur. 
Celle-ci est, dans les cas courants, de neuf chiffres. Ainsi : 1+107!°=—1; 
une variation de 107!° ne peut être prise en compte pour une variable 
qui vaut 1. 
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APPLICATIONS 


Filtre numérique 


Comme conséquence directe de notre étude, nous pouvons essayer de 
comparer les courbes de réponse de filtres numériques et de filtres 
passifs, ayant même pulsation naturelle w,, en fonction de la fréquence 
d’échantillonnage w,. Le programme («Filtrage ») est simple. Il permet 
de tracer les courbes en coordonnées logarithmiques. 

De façon générale, le filtre numérique se caractérise par une suite de 
nombreux Zéros de transmission. En basse fréquence, sa courbe de 
réponse «colle » à celle du filtre passif; le domaine où les deux courbes 
coïncident est d'autant plus grand que w, est grand. 


@6/10 


10 &, 100 





o 


0 dB LEE ee 
| 
be 
| 
| 





20 log|T| 





— 20 dB 


&y =, 


A NT Os 2 mo ul 


20 WO=2%PI/DT 


= 96: 





ÿ à 


nr nent arts arte 





2 
T 
(=) 
[ol 
I 


emmener mtartene verte nettes tente ne a ebenee von 


PI/DT 


20 WO 





nent eona nant non nntansenn mate tn net neutron interet on a permettant near rton 





20 WO= PI/(2*DT) 





ennemi nntemnanlanonennnennnnrenmnnennndennnnnenmnnunnnnmennnnnnennnnnnnnnent 


20 WO = PI/(10xDT) 


[On 





, 


A 
[NS 





20 WO = PI/(50*DT) 


De 


SRE PL UE SE enr 


ki 
| 





es 


Programme “Filtrage” 


14 CTS=IE-S 

24 WASPI.-E SAECT 

34 W=4A.-166 

4 HIRES 

Sù GÜSUE AA 

EN DEF FH FEU el /SQRE +6 AMG TE 5 

ré CEF FH GEL SURC 1 +C 0 TANCHEDT/2 CU 
EDIT 2 Et 

En Hetgcé 1-48 5 


FOR Len TO 259 
W=ll£H 
Fa-sa£tLOGé FH FEU: Ge-SA£LOGCFH Gé 


= 
min 


SR SE EE mm 2 | 


3 rs 
ir 


IF F>194 THEN 154 
CURSETU, F+45, 1 

IF G > 19% THEN 174 
CURSETU, G+5, 1 

MEXT L 

STOP 

PEN AXES 

CURSETA. 4, 1:DRAM 239,0, 1 
FÜR Mel TO € 

234 CURSET H#48-1,9,1:D7 
54 CURSETE. 4, 1:DRAMA, 19 
: DRAM, 193. 1 

258 FOR Mel TO S:CURSETO. 24H45. 1: DRAUS 
«4, L'MEXTH 

264 PATTERN 45: CURSET4RXLOGC Z#P1/€ DT#HA 
»+89, 0, 1 :DRAND, 199.1 

274 CURSETE®, O1: DRAM, 199,1 

284 FETURH 


JM LA Re is Foi Si = DO Si 1 


CE 


St O0 1 


3 PQ res Si à 
ee 
ax 


14, 53 L'HEXTH 
1: ORALE, 0,1 


AL 
ss 


Charge d’un condensateur 
L’équation du système a été donnée dans l'introduction du chapitre : 


RC O4 (H=e(t) 
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Posons e(t)=— constante (cas de l’échelon de tension), s(0) = S0. 

Le programme s’en déduit facilement (voir programme «Circuirc»). 
Quelques remarques cependant. Le raisonnement précédent ne s’appli- 
que pas à l’échelon, qui est un signal non périodique, et donc pour lequel 
le spectre de Fourier contient toutes les fréquences. Une étude graphique 
montre alors que, dans ce cas : 


DT dt 
— soit —<— 
RC " 7 10 


semble convenir. 

Par ailleurs, des valeurs trop grandes pour dt donnent des solutions 
plus grandes qu’en réalité. Cela vient du fait que la dérivée initiale de la 
fonction est prise en compte plus longtemps. Des valeurs trop petites 
pour dt ont pour conséquence une erreur d’arrondi. Un sous-programme 
(300) donne la solution exacte pour le cas où s(0)= 0. 


Remarque 


La solution théorique exacte est : 


—t Const f‘ (u—t) 
s(D=s(0)exp (=) + RC exp ( RC ) du 
0 
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DT/RC—1/100 





s(0)=0 
e(t)=E constante 


Les 2 courbes (solution numérique) et (solution exacte) sont superposées. 


3 courbes superposées 
DT 1 


RC 2 
DT 1 


7 RC 10 


— solution exacte 
s(0)=0 
e(t)=E= constante 





oser serons eg sg eseestnes esse etes 
es RENE e ARE ne en 






eee aan ns 





tentant Prog rmntr ta mmemone romans 0épmrtantané etant cheat nées 66Aane romans ete tement ere 





DT_ 1 
RC 10° 


2 courbes, l’une pour s(0)= — 10, l'autre pour S(0)= +7; on a e(t)=10. 


Programme ‘’Circuirc"” 


INFUT'OT en fraction de RC:":OT 
INPUT" Tension & initiale:"is 
INFUT'AmPlitude de l'echelonitiE 
HIRES 

GOSUE 244 

FOR UsG TO 164 STEF € 4B£CT 
D£stE-S MECT :5=5S4+0S 

 CURSET41+1, 148-SGK4S/E, 1 


CURSET41+1, 148-S4#E-/HESCE 3, 1 
NEXTLI 
STOP 


FEM HXES 
CURSET48, 164, 1:0DRAU 180,64, 1:CURMOVE 
3, 0: CHARESZ, 4,1 


Nm e be if DEN BR D 
ER EE ET RNA EL 


RARE ELN EUX EN 


ul 


= #99: 


4, M: CHARS 4, 4, 1 

34 FORHSITO d':'CURSETA4G+4@#H, 14, 1: DEAN 
4, L'HEATH 

°240 FOR Hez1 TO &: CURSET 36, 2G+C4#M, 1:D 
RAS. 6, 1:MEXTH 

ESA CURSETE. 6,1: CRAUG, 199,1: CRAUZSS, 6,1 
:DRAUG, -199, 1: DRAU-233. 4.1 

264 FETUEH 

244 FÜR L=4 TO 164 

CLRSET dll, € TE PE LL dé EE SE re 1 A0 


24 CURSET4N, 189, 1:DRAU @.-164,1:CUFMO 


CEA 
EE 


24 HEATI 
34 KETURH 


GC ee 


Autre exemple 


La même méthode peut servir à résoudre d’autres problèmes. 
Considérons le cas de la recherche de lignes équipotentielles, c’est-à-dire 
de lignes le long desquelles une fonction donnée F(X, Y) garde une 


valeur constante C. Le programme proposé montre l'exemple (ligne 600) 
de la fonction : 


F(X, Y)=X'—V'+3xXY 


Il demande d’abord les valeurs extrêmes de X et de Y, ainsi que 
l'intervalle de calcul H (soit un At), qui est normalisé par la suite. 
Ensuite, il est demandé une valeur particulière de C, ainsi qu’une abscisse 
de départ A (POKEZ28,0 permet d'écrire dans les trois lignes de texte, en 
bas de l'écran haute résolution). 

La méthode de Newton cherche alors un point initial de coordonnées A 
et Y (A), s’il existe. La courbe est décrite par les méthodes précédentes; 
elle est interrompue dès la détection d’une tangente verticale. La courbe 
est enfin reprise du même point, dans l’autre sens. 


=3%. 2 


FENM ### 
IHFUT 
IHPUT 
IHFUT 
IHPFUT 

IF #E<=z4 


Pl 


HA 


=, 
2 


Li Pire 


es 


=, 


La 
La 


& 1 


Tuuxx: 


Li 


ri mot Li LE OO Ji Cn 


=, 
ee, 


IHPUT In 
CIM “183, 
FOR Kzd TÙ 
YÉéKksiEin 
HEËT K 
HIFES 
POÉE 5: 
ï lignes 
IHFUIT " 
IHFUT " 
IF He 
Mar 
FÜR K = 1 
= À “= 
GOSUE 544 
HET K 
IF .J=f 


ARENA) 
# 


us 
EN 


pe 
Cx | 


—, 
ES 


2 
Lili D CIN int 
Et uit ER] Ex Ü 
iu 
7% ee 2 { 


Et A Ti 
ÉPREANE AE LS 


ii E 


TH 


PU Pie pe he 


Éarai = Ci 
14 GOSUE ESA 
CA IF JA TH 
éasui æ oc": 


K=3 
IF k 
=: 
GÜSLE 
GOSUE 


FX 


D Pat Pi Hi Di 


LE TT Te ES Re di dt 


ET 


un 






AE] 


ez 


— 





FÈR: 
Pa + Pre Pi 


Fri 


4 


= 


ii 1 


:GOEU 


£ IF PENER € 
gi GOSUE Fa 
318 IF x = 0 


SA THEH 


GOT LE 


a 
Ti 


IRC FE 


TENTIEL ##%# 
min Mie 
la Mi 
min MiTh 
Max "ivTE 


OR =YA THEH FIHG: 


*Eé GOT 


tervalle de calcul = "5H 
LE 1 © 
15 


:REN Fermet d'ecrire dan: 
on bas de l'ecran Maute Kes 


an 


ÉHIQN SE Co se Mi 
ee initiale à = "A 


Ê “A THEH # mi 
GOT 154 
DCR -STER ST 
VA + CTE-TAIEES 


FRINT " 


151 


EH FPEIHT " Fas de 
QUOTE 1e 


solution 


EH FEIHT " Fas 
GOTO 124 


dé solution 


1 THEH 34f 


Le Ke 3 


GUSUE 634 


PEAPYAE SORCFE 





E : 


chk4 td: COTO € 


OR #X =260 OR x = à OR 
=-H + Hsk+i: GOTO 835 


La 


ee 
CURE EN] 


ES Li dif 
EURE 
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EN Le ti Pot te ve 
LE 
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EN ELN mm ARRET EL 


CA ER ET ER CA dt Le ai a 


Fat ne EE LG Ts 
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Pr Pr Ti ER CN ER ER 
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LE UE 
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es 
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EX 
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= 
er 
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ee, 


LD Ti in ie ts Mi 
Fi 


m 


Par a 


ÉRARSAE RUE LEE 


Ce cr A tt Ce 


ii CA HR ti 


Jsë: F 
fe Ke 


HEATH 
GÜTO 1€ 


FE H 
GOSUS 
IF FY 
“ Î æ “ 
IF ES 


LTUEH 


ER 


“el 
HET 
ÉETUEH 
RER ÆX 


F = XEx 


FRETUEH 
Mate : 
#=HH-He 


FETURH 
THE 


QE ke 4 T0 
ELG UE SE à 


ns! 


Fi 

= 1 TO 164 

SAM: GOSUE 
A THEH KE 

— F<FY 

A a PU 


# FOHCTIOHS 
HR HUE 


Hed+He {UE 
144: GOSUE 


eee 1 BE 


14 
sf: Li ESkE+1 et 


San 
TURH 


THEH Je: YCEÉ e#l 


F, Fr, Fu ŒEY 
+ SERA T TC 


GOSUE GGg:FHeF 


SAGE: Fac FH-E CAUSE H 


SOEUR CES: FHer 


YeYH-H/ 180: GOSUE AD: Fac FH-F ESC 


“H 


} RETURN 


REN #E 
LEE 





£ TRACE &X4X 
me Ce SE TAELS CE ne mt Et 
TH 
# 16 THEN a = LE 
CA THEN VX = 
+ CA THEM 6% = CH 





FETUEH 


su THEH fre 


Î 


ua True Es 


RE ÆXX Valeurs 1 


ksg: J= 


et 
pe 


J=.+ 1 : 
Hak4+ 1: 
IF CE 


Hi: IF VCGe 
LI EX l'une 
IF K=11 TH 
12TE OR YCE 


FÜR Leg TO K-1 
FF CCR AT TELE 


bte PES. JD UTHER 
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EH FRETUEH 
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3. MOUVEMENT PONCTUEL NEWTONIEN 








Un mobile ponctuel, soumis à une force centrale centripète d'intensité 
inversement proportionnelle au carré de la distance, décrit un mouvement 


newtonien. 





T=OM 
u : vecteur unitaire 


Intéressons-nous au cas d’un satellite de la terre, que nous assimilerons 
à un objet ponctuel. On peut écrire : 


MmMy — 





F=ma=-G— 
T2 


Cette équation conduit à un mouvement plan. 
En appelant G(—=6,67 X 10 !! MKSA) la constante de gravitation 


universelle, on aboutit à : 


Gmry — 
à 





a = — 


T 
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Il s’agit d’un système d'équations différentielles, que l’on peut écrire : 








d? x x 

= OGM 
dt CC + 7) 
d? 

RE ue 
dt (x + y) 


La résolution de ce système dépend de deux conditions initiales (à 
t=0), on se donne en général T, (distance au centre O) et la vitesse de 
satellisation vi (module et angle de départ). 

Il y a (au moins) deux méthodes pour résoudre ce problème : 


e Une résolution directe du système d’équations, en étudiant l’évolution 
du système pendant des incréments de temps dt. Cette méthode a 
l'avantage d'utiliser ainsi le temps comme paramètre, ce qui permet de 
suivre le satellite en temps réel. 


e Une solution analytique du phénomène, c’est-à-dire le tracé direct de 
la trajectoire. Le temps n'intervient alors pas. 


PREMIÈRE MÉTHODE 


Elle est classique et a déjà été étudiée dans cet ouvrage. Nous posons 
K—Gm.,. L'intervalle de temps choisi est DT— 10 secondes, ce qui 
paraît raisonnable. Les conditions initiales sont : l'altitude d’éjection H, 
l'angle d'éjection A et la vitesse d’éjection VO. 

Le programme évalue dans un premier temps (à titre indicatif pour le 
lecteur) la plus petite vitesse permettant la satellisation dans le cas : 

T 
LT 
Un facteur d'échelle F est calculé, prenant le rayon de la terre R, 


comme référence si G—1. 50 points son l'écran représentent R,. Si 
G—0,5, cette échelle est divisée par deux. 


Les grandeurs initiales en coordonnées cartésiennes sont ensuite 
évaluées. 
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Le système d'équations différentielles est enfin intégré. Trois exemples 
de courbes obtenues par ce programme sont fournis; on notera que, dans 
le troisième cas, le satellite s'écrase sur la terre! 


19 K=9.S#C6.8"2#k1E12 
24 INPUT'ALTITUDE C'EJECTIOH Cen m ="; 


H 

34 DCT=1g:Ha=Eé, SEE 

4 MMSSORE ZEKKRA/E C RAM EC LERAHH 5 1 

54 FRIHTYM 

6 IHPUT'UITESSE C'EJECTION Cen ms =" 
if 

F4 IHPUT"ANGLE C'EJECTION £en dedres = 
OL :H 


24 DCE D'ECHELLE ='";:G 

24 FeG#SA/'é. SEE 

LA FzFA+H : fes MALCOSE AKLEPTI. REA 1: Me UTE 
STIHEMEZEPTI-GÉéG ee RE4H VE 

114 HIRES 

124 CURSET1I20., 144, 1:CIRCLESAEC, 1 

154 CURSETA. 6, 1: De. 1:DFALUG, 199,1 
: CRAU-ESS, 4, 1: GRAUR, -199, 1 

LA RESRTS : He ERA: AVE -HOEN RS 

154 Dre HOT : EMA EÆCT 

UP ET 2 Me I ECT 

CURSET LEA EGHER 1, LACGI-CNCER 1, 1 
Fe SRE GENE UN 
TISTI+DT:TEIHTÉTI SEGA: FRIHTT: MH, 


DE Jin 
+ 


Loeb pe po 
Br ji 
ÉNSCRE OC EE EN EL 


FREIHTT1 

GOT 4 

LFRINT'FACTEUR C'ECHELLE = "56 
LFRIHT'ALTITUCE GE CEFHETU: Hi" mm" 
LPRINT'YUITESSE D'EJECTION":4V@:" ms 


Pa Pet Fo Pa Po 
5 Pa dv 15 


ENSES 
RE ri 1 


if 
in 


LÉRIHT'ANGLE C'EUJECTIOHNM: Ai" degres 


264 LPRINT'YITESSE MINIMUM =";9M;" me" 
274 LPRINT'PERIODE CE REVOLUTION =":T1: 


2 STOP 
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SECONDE MÉTHODE 


On démontre que la courbe décrite a pour équation : 


1 _Gmr 


—+E cos 0+B sin 0 
T € 





avec : 


C=7T, vo sin & 


l Gmr 
Es DES ESC ENS 
To To Vo SIN 
Cotg @, 
To 


Les conditions initiales à introduire sont les mêmes que précédem- 
ment: Ty, Vo, À. S'agissant de la solution exacte du problème, la 
précision du résultat est liée au choix de l’incrément de l'angle 6; 
l'équation étant en effet du type polaire implicite T = f(0). Les trois 
mêmes exemples ont été choisis, à fin de comparaison. 


14 K=9,8#(6,8"20#1E12 

2 INPUT'ALTITUCE C'EJECTIOH éen m =": 
H 

39 DT=18:FG=6,8E6 

49 VYM=SAQRC 2HKKRO EC RGEH IEC DEREHH 3 

24 PRIHT#YHM 

64 IHPUT'YITESSE C'EJECTION Len m/s =" 


: a 
74 INPUT'ANGLE D'EJECTIOH Cen degres De 
":A:A=20-A 


89 INPUT"FACTEUR C'ECHELLE =":G 
QG F=<G#54/6,8E6 

144 R=RG+H : DSREVBESTNE AXTEÆPT 360 1: EC C1 
SRI K SC DEC 30: Bz-1/CRETANCAKEXFI/3606 3 

114 HIRES 

124 CURSET124, 140, 1:CIRCLESAKG. 1 
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134 CURSETA. 4,1: DRAUZ3S, 4, 1: DRAUG.199,1 
:CRAU-239, 04,1 :DRAUG.-194.1 

144 T=g:C2=0#0 

154 T=T+,.45 

164 Relié K/C2 146 EKCOSCT 34e EXSTHOT 0 
174 H=FEREACOSET 5: YEFERESINHET 

180 IF K>119 OR #€-119 OR “2:99 OR v<-99 
THEH 244 

194 CURSET 124@+%,1800-"%, 1: GOTO1SA 

ZA INPUT'Vaulez vous continuer ":A$ 
214 IF A$ = "O0" THEW GOTO 15 ELSE STOF 





Programme 1 


Facteur d’échelle =] 

Altitude de départ 1000000 m 
Vitesse d’éjection = 8 000 m/s 

Angle d’éjection 0 degrés 

Vitesse minimale —= 7356,43964 m/s 
Période de révolution = 7 540 s 
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Facteur d’échelle 
Altitude de départ 
Vitesse d’éjection 
Angle d’éjection 
Vitesse minimale 
Période de révolution = 7 520 s 


Programme 2 


= 1 

1 000 000 m 

= 8 000 m/s 

10 degrés 

— 7356,43964 m/s 








Programme 1 
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Programme 2 
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Programme 1 


Facteur d’échelle = 0,75 

Altitude de départ 1000000 m 
Vitesse d’éjection = 8 000 m/s 

Angle d’éjection 20 degrés 

Vitesse minimale —=7356,43964 m/s 


Période de révolution = 0 s 
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Programme 2 
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4. RÉGRESSION LINÉAIRE 





Considérons un ensemble N de points, de coordonnées X et Y. Le 
problème est de trouver une droite passant entre tous ces points, de telle 
manière que la somme des carrés des distances des points à cette droite 
soit minimale (méthode des moindres carrés). Si tous les points sont 
alignés, il se trouvent évidemment sur cette droite. On dit alors que leur 
corrélation est parfaite. Inversement, dans le cas général, il est possible 
de définir un coefficient de corrélation s permettant de chiffrer dans 
quelle mesure ces points sont alignés. On démontre les relations 
suivantes : 





14 INPUT'Hombre de couples de Foints:"i 
WN 

24 Sk=û : 5 Y=A : 220: 220: Ll=@ 

34 FOR J=1 TO H 

48 PRINT. : INPUT ia 

45 LFRINTI, 4." 
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54 £ 
EM & 
74 & 
84 Y 
ag 1] 
149 
114 
124 
134 
144 
145 
w + 
154 
164 
174 
189 
set: 
134 
135 
t de 
298 
214 
31224 
2ef 
238 
23% 
M AYX 
244 


En posant : 


RARES 
Y=S +" 
LRALHAKE 
L=E+YE" 
UE" 
HEXT J 
sé NELI-SHKS NY IE HEAZ-SHAESE 
Es SY-A4SEA I/M 
PRINT : PRINT :PRINT 
PRINT'YG x = MA" x + UE 
LPRINT:LPRINT:LERIHT" Yémen = "Hi" 
Ve UC GREEN 0 
W=v2e6 SES NS 
S=V/1 
PRINT'Le coefficient de correlation 


PRINT" sg = "SR S 
LPRIHT:LPRIHT:LFRINT" Le coefficien 
correlation est :":SORCS 0 
PRIMT : PRINT :PRIHT 
FRINT'Voulez-vous estimer Un Ÿ Conf 
nt #7": GETAS$: IFAS="MUTHEN STOP 
INPUT "sem"; x 

PRINT", = "AXX+E 
LPRINT:LPRINT:LPRINT" Y est egal 4 
+Bs"pour x = "i4 

GOTQ 219 


N 


SX = somme des X=X X; 
1 


N 
SY —somme des Y=Z Y; 
1 


N 
X2 = somme des carrés des X=X X° 


1 


=, AT.= 


N 


Y2—= somme des carrés des Y=Y Y? 
1 


N 
U —=somme des produits XY=X X,Y, 
1 


on obtient : 


__N:U—(SX-SY) 


1 
= : B=— (SY—A:SX) 
N-X2— (SX) N 





Sx-SY 
1) 


s? = 


Ris 





EXEMPLE 


Considérons l’ensemble des températures moyennes relevées à Mar- 
seille en septembre 1984. La coordonnée X représentera le quantième du 
jour concerné, la coordonnée Y, la température moyenne de ce même 
jour. La loi obtenue (voir le tableau) est Y = —0,16 X+23,6. Le 
coefficient de corrélation, cependant, n’est pas très bon : 0,66. (Cela est 
dû au fait que le temps a des «caprices » en cette saison de l’année.) On 
peut déduire de tout cela que : 


e la température moyenne de septembre a été : 21,13°C (X=15); 


e la température moyenne de septembre a décru de 23,43 °C (X=1) à 
18,66 °C (X = 30), soit une chute de 4,77°C au cours du mois; 


e la température moyenne d'octobre sera (X—45): 16,19°C, cette 
prévision étant toutefois faite sous toutes réserves! (En fait, elle a été 
de 18,80 °C.) 
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L 1 24 
3 3 24,5 
4 4 25 
5 pe Le 
6 6 13.5 
7 7 19.5 
8 8 24,5 
S 9 21.5 
14 19 22 
11 11 22 
12 12 22 
13 13 21.5 
14 14 22 
15 15 23,5 
16 16 22 
17 17 23 
18 15 24.5 
13 19 18.5 
24 20 21.5 
21 21 21 
22 22 19,5 
22 23 15 
£4 £4 17 
25 25 17 
26 26 18 
27 27 21.5 
29 LE za 
24 24 22 
34 34 19.5 


Veux) = —.164516129 %X + 253.6 


3, 425482 pour 


14.6645161 Pour % 
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5. Quelques fonctions importantes en analyse 
6. Dérivée numérique et chaînage de programmes 


7. Intégration 





5. QUELQUES FONCTIONS IMPORTANTES 
EN ANALYSE 





INITIATION AU PROBLÈME 


Un signal modulé en fréquence est de la forme : 


u =U sin(Q, 1 + g, + x sin wt) 


x indice de modulation. 
Et si @=0 : 
u = Ufsin {,, t cos(x sin wt)+cos (, t sin(x sin wt)}. 
On peut développer les termes : 
sin (x sin wt), cos(x sin wt) 


en une série de Fourier dont les coefficients sont les fonctions de Bessel 
de première espèce : 


nn | 
Jo (x), Ji Cx), D Le F(n+k+1)k! () 


où l(x) désigne la fonction gamma : 


F(n+k+1)=(n+k)! 


On trouve finalement : 


u=U{sin Q, 1(J,(x)+2]J,(x) cos 2wt+21J,(x) cos 4 wt +...) 
+ cos (, 1(2J, (x) sin wt+27J,(x) sin 3 wt+...)} 


=U(J, (x) sin Q, 1 + S J, (x) {sin(Q,+ now) t+(—1)" sin({Q,— no) t}) 


0 es 


et l’on voit aisément une infinité d'ondes latérales dont les amplitudes 
sont liées aux fonctions de Bessel de l’ordre n et de l'indice de 
modulation x. 


3 x 
Les ondes latérales de rangs supérieurs à Eu sont négligeables. 


La Figure 5.1 montre quelques fonctions de Bessel. 





Figure 5.1 a : Fonctions de Bessel. 





Figure 5.1 b : Fonctions de Bessel. 


69e 
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Figure 5.3 : Fonction d'erreur. 
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MÉTHODE 


On propose donc plusieurs programmes qui servent à calculer les 
fonctions de Bessel et la fonction gamma. 
Le premier programme est basé sur le développement en série : 


x\" = 
_ 2 b 
(:) k=0 . 


J(x, n)= 

F(n+1) 

— b,_, x? 
bo=1, b=—"——, k=1, 2, 3, … 
à KE (4k(n+k)) 


Il faut introduire le paramètre n, la valeur F(n+1) si n n’est pas 
un nombre entier, et la valeur de la variable x 0O<x<7+}n|. La 
sommation est tronquée selon la précision voulue (sept chiffres). Le 
deuxième programme donne aussi les fonctions de Bessel de seconde 
espèce en se servant du développement asymptotique : 


J, x) + jY, x) = 





1 
r(m+n+-) 
= 4) ei-7/2n—7/4) Les J ENS RO LE 
TX m=0 (2 x) mr (m-n+s) 
2 
1 
SH r(k+n+s) 
J 1 
+) lel<1 n>—> 


Pour x=7+|n|, sept chiffres significatifs sont garantis. 

Quand 2 n est un nombre impair, la série ne possède qu’un nombre fini 
de termes. Dans les autres cas, on arrête la sommation dès que les termes 
recommencent à grandir. 


20562 


Le troisième programme concerne la fonction gamma. Nous rappelons 
sa définition : 


rG= | e-tt—ldt, u>0 
(0 


C'est pourquoi une intégration par parties donne : 
F(u+1)=ur(u), F(1)=1 


et on déduit : F(n+1)=n! pour tout n entier = 1, 2, 3, … 
En se servant de la relation : 


F(z+m+1) 


D ER | EE PE 


on voit que la fonction possède un prolongement analytique dans le 
domaine : 


€C/{0, —1, —2, —3,..} 
On a donc F(n)=co pour n=0, —1, —2, —3,… 
Pour évaluer F(u), on utilise son développement asymptotique : 


1 1 1 1 


1 
L'(u)& V2 ( =)  — 
(u) resp {(u DOM RES ET cou "lo LœUe 


+ —.}(u26) 


On peut attendre sept à huit chiffres significatifs (Figure 5.2). 
Un quatrième programme étudie la fonction d’erreur : 


7 x 
o(= VE | e-"”2dt 
0 


LG: 


Dans le but de calculer la distribution normale, on utilise sa densité de 
probabilité : 





La fonction d'erreur est donc étroitement liée à la loi normale (voir la 
Figure 5.3). Dans l'intervalle [0,3] on se sert de son développement en 
série : 


4 6 


EX —x2/2 [1 w “ + Lo 
2e ser + +) 








pour x>3, de sa fraction continue : 


D(H=1—— e-*° (He + ) 
V7 X X X X X ù 


On tronque les sommations après 23 (resp. 16) termes pour garantir 
huit chiffres significatifs (!). 


ORGANIGRAMME 


Comme les méthodes des quatre programmes proposés se ressemblent, 
on ne donne qu’un seul organigramme concernant la fonction gamma. 


1. An a priori Estimate for the Truncation Error of a Continued Fraction Expansion to 
the Gaussian Error Function, G. Boesse, Computing 29, 1982, p. 135-152. 
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Programme : Fonction gamma 


calcul de F(U1) par 
son développement asymptotique 





LES PROGRAMMES 


EEN ##4 Fonction d'erreur ### 
REN # IE SORE SCPI AE ARE ETES GK 


& is 


fr 
a it 


Pû 
[ER] 
it 


INPUT "Waleur de la variable x = " 


Pa PO F9 I Po POUR 9 Xe Phi 


234 IF 4 > 3 THEN 2388 

24m S=1:B=1 

264 FOR M=3 TO 45 STEP 2 

264 B = B#A#k-tt: SeS+E 

ra HEXT 

BEN S = SHMKDORCZ/PI MHEHPCK#K/2 
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0 RFI X Po 
Let ai da te Pi 


at Pire 0 = NS 


ru f 
Fi 
L 


3 

LORERR 

cn 

St 1 


RARE 


248 
2485 
2416 
2415 
P42t 
2456 
2435 
Ed 
2445 
2454 


FEIHT "ICSORE AS /PI TÉKEXPC-HTE 2 , 4," 
= # 15 4 EHC 

Ex 

FOR M=15 TO 1 STEF -1 

Se + M/S 

HEXT 

Ga 1 SRE PIE EMPE REA IC ES 
GOTO 2294 


FENM XX Fonctions de BESSEL ### 
REM ÆKX Jun. Programme 1 %## 
IHEUT "Parametre nm = "M 

IF HESINTÉHS OR HÉG THEH Zd445 

Gel: FOR KE=û TO H-1 

GG K +1 1 

HET 

GOTO 45 

INPUT "Garmaé ati + "ic 

PEN EE Woir Frogramme fonction 4 


arma XX 


24 
1,4 


def 
Las 
INT" 
s474 
car 
2494 
20 


DU 


DM 
mt 


in 


M 


PO + 
io ii 


Fr Pi 


HR 


IHPUT "'aleur de la variable x = " 


RENM ÆXX MéxeT+abs( ns #XX 

IF Céée@i Of Che r+HBS Mis THEN FH 
ve 7 4 abefnii voit Programme & " 
BEM ### Hlaorithme ### 

n=i: Mel: JE 

E RE d'AITKE M4 1 

3 J + B: tisM+i 

IF AESCE :1E-4 THEH 2444 

JE Ke 253" MAG 

PRINT MJC"; MMM" = "3: END 


FE #4X# Fonctions de BESSEL ### 
REIT ŒHAX sont. Wéxini Frogr, ZT ## 
EN ZX DeveloPhement asamPptotique 
#XX 

PEN #Y# Precision: 7 chiffres ### 
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£ASA INPUT Parametre on = "5H 

2éd4n INPUT "Waleur de la variable x = " 
4 

2ESA IF € PAHBECNS THEN PRINT x trop 
Petit ": IF #<=4 THEN FING: GOTO Zédf 
PGA REM ##4 Alaorithme #4% 

ETAA B=1i: M=l: Rzi: I=û 

er14 Cab: BSEXC -, SM IT KE SM EME 
#2 


2724 IF MSIHTOM-ZUEZ THEH FR = KR4#f-13%cn 
#81#E: GOTO 2748 

2739 I = I - C-12"INTM/ZI4E 

2744 IF CE=<G3 ORCCABSCB 32 ABSCC 13 AMD € 
M> AES CH553 THEH 2764 

cran M=M+1: GITO #18 

2ré4 = SORCRER + IXI5 

270 F = ATHCI-R5 

2754 IF K< 4 THEN PF = FI+P 

EPA4 J = SURESCPIEE 2 IAAXCOSE K-H#PI/2-P 
1-4+F 53 

LS = SURCZ/C PTE 3 IEAASIMCX-HXFI/2-P 
I#4+P 2 

2819 PRIMT MJC": Mi Min 
ZEEG FRINT MONS Hi N Min 
2634 END 


Hi 


HE 


4 
Hot 
4 


ca4a REM ### Fonction GAMMA ### 

2914 IHPUT "aleur de la variable u = " 
il 

2324 IF = INT ANG U£et THEH PRINT 
"CANMMAS "Un" & infini ": END 

2444 FENM #kk Alaorithme #X4## 

2459 LET Fzi:LET Uiall 


2364 IF Lé=S THEH PeP4U: LET Usl+1: GOT 
2364 

cAr4 LET Wae1fUI#U 

cars LET YaiéC, PSÆV-A VEUT, St 1 KM 1 1 
#-12+1 

2980 LET GESCRC CÆPI THEXPE CUS, 5 MELHOUT — 
LEE" 22 
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5 IF Ui = IHTCU1) THEN G = INTCG+.5) 
9 PRINT"'SAMMAC "UT; ") = "iG 
= END 


EXEMPLES NUMÉRIQUES 


A l’aide des quatre programmes nous avons dressé un petit tableau : 


J, © 


nyx ! 2 3 4 5 6 7 8 9 F(n+1) 











































A|arr E -H.4  -ÿ 4.15 H,1r -1.43 1 
14.44 un. £ -H, [e F1 ï 23 WE 1 
2 | .i1 E £ E 4, Ûs [a] -.11 H.1 ë 
BU GUZ 6. CR f, Ë A, A, -H,23 -1,1: ë 
4[ x a. 5, &. (s] 5 A. -H.11 -5 =4 
514 [5 H,ïd4 E E ñ st H.12 — 129 
CS A5 ä 4,81 4,895 4. a, La) H,ï4 réÿ 
FE a ä BU W, H.15 W,Z ä 141 
ER A ä ä ä 5, H,86 4.13 4. ï 
3156 ï (s] ñ LS) A.HZ M.HS 4 [5 





GAMMAS 16 
GAMMRL 5 2 
SANTA 





ction d'erreur 
G, .E 
9, .4 3 
D: .6 2 
D, .5 5 


AE! 





1 
+= = 
it 










dB 142 3 
nb, 14 à 
n, 1,6 
5, 1.8 à 





4 ES 
M, où = 
A, 4 = 
do = 
TE 











6. DÉRIVÉE NUMÉRIQUE ET CHAÏÎNAGE 
DE PROGRAMMES 








INITIATION AU PROBLÈME 


Nous avons déjà introduit les fonctions de Bessel de seconde espèce 
Y,, (x) et proposé un programme pour calculer leurs valeurs si la variable 
x n'est pas de faible valeur. Ce dernier cas sera traité maintenant. 

Notons d’abord la relation : 


J,(x)cos nm —J,_,(x) 


Y,,(x)= n £int(n) 


sin AT 


Nous trouvons donc facilement Y,, (x) à l’aide des fonctions de Bessel 
de première espèce (voir Fonctions de Bessel-programme 1) si n n’est 
pas entier. Sinon, on applique la règle de l'Hôpital et on obtient : 


1 9] — 1)" 
vo = dr) LD" a, 


s - sm=0, +1, +2 
T Ôôn |[n=m T on |n=—-m 


Pour ne pas trop compliquer la discussion, nous posons m = 0 : 


y. 0 = 2 0) 
X)= = 
T ôn |n=0 


Notre premier programme permet alors le calcul de la dérivée 
numérique d’une fonction et, par un chaînage de ce programme aux 
programmes «Fonctions de Bessel» et «Fonction Gamma», on obtient 
Yo (x). 


62° 2 


9 — y 


Remarques 


1. Rappelons la notation de Landau; on dit : 


f(h)=0O(g(h)) 
si le rapport Î reste borné par une constante au voisinage de h—0. 
8 


2. L'ordre O de la méthode correspond au nombre d'itérations. 


MÉTHODE 


Expliquons d’abord le programme « Dérivée numérique d’après Rom- 
berg ». 
Nous notons les développements de Taylor de f(x+h)et f(x—h): 


FE = FOR CO Hp fo +È À F0) )+O(h°) 


D'où l’on obtient les différences : 





pbs CERTES) h, = f keto 
2h, DEA 
= fo + © + o (5) 


et, par une première extrapolation : 


D(k+1)—D(Kk) 


4 
+ OH) 


f(X)=D(Kk+1)+ 


En réitérant ce procédé, nous trouvons le schéma suivant : 


f(a+h)—f(x—h) 


PES 2h 


263 = 


De D). DQ)+2 2-20 
h h 
flx+=)- (>| 
4 4 D(3)—D(2 
MP Eu dt DG) <— D +220) 
2 
D(2)—D(1) 


D(1) <— D(2)+ _ 


= f'(x)+0(h°) 


etc. (selon l’ordre de la méthode). 
Nous continuons jusqu’à ce que deux approchées successives D (1) ne 
diffèrent plus (précision PR). 


UTILISATION DES PROGRAMMES 


a. Dérivée numérique 


Il faut introduire : 


la valeur de la variable x : X =? 
l'ordre de la méthode o : O =? 
la précision PR : PR =? 


l'intervalle de calcul À : H=? 


Pour réduire l'influence des erreurs d’arrondi, on choisit : 


ne 
IA 
T 
IA 
D 
a 

[ee] 
IA 
© 

IA 
J 


OU : 


nl 
A 
© 
© 

g 
IA 
Ô 

IA 
© 
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b. Programme «Chaînage » 


Se basant sur la définition : 


v, 2 1m 2%) 


T n—+0 n 
on introduit la valeur de la variable x et on calcule : 


dJ, (x n) 


=:J 0 
ne V0 n (% 0) 










Dérivée 


GOSUB 2400 
Calcul : 


J(xh,),I(x —h,) 


Programme Dérivée 


Sous-programme 1 Bessel 





GOSUB 2910 








Calcul : Û 


Sous-programme 2 Gamma 
Th +1),T7(—h+1) 


LES PROGRAMMES 


14 RENM ##Y# Derivee numerique *#%# 

24 KEU #4 d'apres ROMÉERG 

Qn REM #X# Fonction fx L##2 

184 DEF FH FX = ATHC 40 

114 PRINT " féxiu = arctan x " 

199 REM #4% Valeur de 14 variable x #4#%X 
24 IMPUT " f'ésie variable x = "Mix 

214 INPUT " Ordre de 1a methode oz "it 

228 IF <O<é=@5 OR CO<YINHTÉO NY THEN PIHG: 
GOTO 216 


24662 


cs INPUT " Precision FR = ";:PE 

240 INPUT " Intervalle de calcul Hz "5H 
254 IF CH£=G5 OR CHi.Si THEH FIHG: GOTO 
244 

264 Din DéCO+15i 

ce PRINT " Tableau des aPProchess :" 
294 FEM ### Alacrithme *X## 

349 DEF FH DCH5i = CFH FC — FH FE 
NUS EXH à 

314 C1 = FH DH: Di=tc1i 

324 PRINT "K = :",D01 

33A FOR K =1 TO! 

544 H=H-E 

354 PRINT "KÉ=MK 

264 CékK+15 = FH DCH) 

Sara KR=il: RENM ### ExtraPolation ##% 

384 FOR K1 = K TO 1 STEF -1 

299 R=4#R 

44 DK) =DÉK1+10+ CDCKI1+15 — DÉKI Hi 
R-1 

414 PRINT CéK102: 

4ZA MEXT K1 

428 D2 = Di1) 

440 IF HBSD1 - Deiéz ABSCPRI THEN FRIN 
T:PRIHT " Cerivee f'CMiki" sn æ "iD2: END 
454 D1 = Ce 

4a£A PRINT: HEXT K 

478 PRINT " Precision "iPRi" non obtenu 
e,. augmentez © " 

459 PRINT " Derniere aPProchee : " 

434 PRIHT " f'<":4i"5 = "iDg: EMD 


5 REM ##% Chainase +2 

14 FENM ### Derivee imerique ### 

“9 REM ##* d'apres ROMBERG EE 

114 FRIHT "fénis Jéx,n2" 

115 INPUT "Valeur de 14 variable x = "; 
# 

124 IF SKké=4n OR Cke=F4+ABSCHii THEH FAI 
HT! be 8 + abafnis voir Frogramme 2" 
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124 REM #44 Valeur de 14 variable x ### 
214 INPUT " Ordre de la methode oz "if 
234 INPUT " Precision PR = "PR 

£a INPUT " Intervalle de calcul he ":H 
264 CIM CéO+i 

284 FEINT " Tableau des aPProcheezs :" 
2939 REM ### Algorithme #4#* 

294 H=H:GOSUR 2444 :.J1z=J:Hz=-H:GOSUERZ4GA 
214 Déti = Use e#Hi: CisDc 15 

SA PRINT "K = 4 :", 01 

334 FOR K =1i T0 0 

344 HEH-E 


LG PRINT "K= "ik 
255 HeM:GOSUB244û:.11=.1:H=-H: GCOSURE 444 
SEQ Dék+1n = CJIi-U1/C2#H 5 
74 Rel: FENM #£# Esxtralolation ##* 
364 FOR K1 = K TO STEP -1 
394 Fzd£&R 
ann Dékiin =DCK1+1 4 CDOCKÉL+1n — DT IeE 
F1 5 


414 FEIHT DCKi. 
24 HEXT Ki 
438 DZ = ti) 
448 IF ABSCDI — D23<= ABSCPRI THEH PEIH 
PRINT" Derives Jné"i#i"@n = ":02:GOTO 


in 71 


SA 

454 Di = DE 

dé PRINT: HEËT K 

4FA PRINT " Precision "FRS" non chbternu 


. augmentez co" 

488 PRINT " Cerniere aPfProchéee 1" 
AA PRINT "Et imite Mi EKOUP-PIT: END 
sd FEI ### Fonctions de BESSEL #4 
CAB REM XX xt, Programme 1 Æ## 
cit REM 4 in Programme fonction SG 
na AK 

LASA GOSUR 2914 

2A4FA RENM ### Alacrithme #*% 

ARR Bel: Mel: SE 

2494 À = EH H ce D'EMTEC HAT 

LEA JU = I + FE: Mets 

ELA IF HÉSCBILIE-S THEH 249 


il 


e 
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JE NEC SAINS 

FETUREH 

FENM ### Fonction Gamma HW+1i in ##% 
RETSESI 

LET Fel:LET =] 

IF L£SS THEH P<zP£U: Uallti: COTORSE 


ai fai Pi EN fo PO 


OR IDUENA 
Di ee 1 Lg Po 


EN OU RAOCE LOUE OURS | 


Der 


2374 LET Y=i<CUÆU 5 

ESEG LET VecCc PAM EME. Si DEV 1 
HU 1241 

CRÉES GESQRC SEPT IHESPE CUS ÉLNOU LE" 


EU - 
iû 
LD 
Si 


fi 


FRETLURH 


EXEMPLES NUMÉRIQUES “DÉRIVÉES"’ 


FX = arctan 

- variable x æ % 1 
de La methode Qi 
Fréecision FF = % 1E-6£ 
Intervalle de caleul he FT ,5 
Tableau des aFFrochess : 
E & H : 15146114 
E = 1 


OUEN SE.s 





ee LT 


EH = ZE 
OASIS TE  . AGASNESISS 
H = = 


Exemple 2 





Lo HS QE mets 1 1 
| variable x = Le. 
Ordre de La methode ; 


it ir 
à 
Eu 


uk ÿ 


Precision FR = % 1E- 
Intervalle de calcul he % ,5 
Tableau des aFFrothéees : 

HE = A : OASIS PRRSE 


F 
lu 
ne 


EH = 5 

ST Perl  .SPFESUSES ,Srraotedé 
.Srraselid SP FISUET Sfr soneé 
SAC ele 

Deriveg F'C E = Sr 35 

valeur exacte: .SPPSSOET 


EXEMPLES NUMÉRIQUES ‘CHAÎNAGES" 


Nous allons calculer quelques valeurs de la fonction de Bessel de 
seconde espèce Y,(x) : 


Fous IC na 














Malgur de La variable x & 1 
e de la méthode ce d 
cision FR & 16-68 
Tuterealle de calcul He SE 
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Me À 


Fénie sine 


Ma leur “4 & | 












ter 


1 LEE X 









Der. 
“ES s 1 L 


a 1 ire hp . 


IE 


— 70 - 





Po 
Lac P EN 





7. INTÉGRATION 





INITIATION AU PROBLÈME 


Partons de l’équation différentielle d’un objet en rotation obéissant à la 
loi de la pesanteur (voir la Figure 7.1): 


mL?7+Lmg sin z=0 


mL? moment d'inertie 
z(t) angle L 
Lmg sin z couple de rotation 


t temps. 
mg 


Figure 7.1. 


Nous trouvons facilement l'intégrale première de l'énergie cinétique en 
intégrant : 


Ne 


= 28 
2=— cos z+C 
L + 


Admettons les conditions initiales suivantes : 
zZ(0)=2z» 2(0)=0, —T<z,<0 
L'objet va donc descendre et il en résulte, pour t de faible valeur : 


_ 2 g(cos z — cos z,) 
_ L 


°2 


72 


71 


FE fé dy 
= VE f = — 
28/., Vcos y—cos Z 


ff ee = V7... | 
LEEDS : 


On ne trouve pas de primitive pour évaluer cette intégrale. De plus, la 
fonction à intégrer présente des singularités à y= + Z4. 

Nous voulons d’abord écarter ce problème supplémentaire et étudier 
une méthode pour calculer l'intégrale définie : 


ef f(x) dx 


d’une fonction régulière (continue et bornée) sur l’intervalle d'intégration 
[a b], — oo <a<b< oo. 
Cela nous donnera l'écart de temps : 


Dt=t(z)—1t(z), 2 <a << —2%Z 


que l’objet met à passer d’une position z, à l’autre. 


Figure 7.2. 
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MÉTHODE 


Nous choisissons la méthode de Romberg et gagnons une première 
valeur approchée à l’aide de la méthode du trapèze : 


14=1@+ 0 h h=b—a 


Ensuite nous partageons l’intervalle [a, b] d’une façon régulière : 


h 


h <«— - 


et nous trouvons de même : 





1@)=(fta)+25 (2) +50) 
En combinant : 


1@)=ECT) 


1(1) <— I1(2)+ RE 


on obtient la formule de Simpson bien connue et souvent appliquée : 


le 2 (rc +45 (TE + 560) 


Procédé général : 


(b— a) 


Arrivés à h = TL »k=1,2,..,n(n : ordre de la méthode), nous 





posons : 


h 241 
LH DE (F(a+ F4 E f(a+ih)) 


2k-1—1 


= (+ h ZE fa+Gi+Dh) 


= 7% = 


et ensuite pour k1=k, k—1, …., 1: 


((k1+1)—I(k1)) 


PER NEED ae p 


1(1) nous livre une meilleure approximation de l'intégrale I. L’erreur 
I—1(1) aura la forme : 


(b—a} "+ B;,,+2 


1—1(4)=( 1 ————— 
= ) 2+0D  (2n+2)! 


f2r+2(C) 


(B, nombre de Bernoulli, Ce[a, b].) 


ORGANIGRAMME 


Intégration d’après Romberg 


Entrées Fonction f(x). 

100-285 Intervalle d'intégration [a, b], L= b— a. 
Ordre n de la méthode. 
Précision relative P. 


L b 
Initialisation  1(1) <— QE FE) 
300-320 I <—I(), M ; 


L 
Re RER 


a. Calcul I(k+1): 


k1=1,3,5,..,M—1 

SR SRE RLR) 
I(k) 

KÉÇUETTSAR 
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b. Calcul d’une meilleure approchée I (1) : 
k1=k, k—1,...,1 
R «— 4R 


A(k1+1)—I1(k1)) 
(R—1) 


Ta ICT) E=-T> OUI —+ FIN 


NON 


c. Précision relative P obtenue? 


[12—11|< Px/12[? 


NON M «— 2M,I1 «— 12 


=> OUI ——+ précision 
non obtenue 


I(k1)=1(k1+1)+ 


OUI —+ I1=12 





FIN 





NON | 
k «— k+1 FIN 
Variables 
A, B Limites (inférieure et supérieure) de l'intervalle. 
L L=B—A. 
N Ordre de la méthode. 
P Précision relative. 
R Variable d’incrémentation. 
K,K1 Variables de boucles. 
I(K), I(K1), 
11,12 Valeurs approchées de l'intégrale I. 
M, H Subdivision de l'intervalle, M=2#, H == 
S Somme. 


= 752 


PROGRAMME 


19 FEM ### Integration #4# 
24 FEM ### d'apres ROMBERG ### 
GA FEM.XK#E Fonction FCI ÆKX 


144 
114 
199 
aa 
14 
ce 
239 
z4@ 
234 


GEF FH FL4D = #77 

PRINT Mféxi = KP" 

REM ### Intervalle Ca. HT *## 
INPUT "Limite inférieure 4 
INFUT "Limite superieure b 
IF B<£=A THEN PIHG: GOTO 244 
FRINT MECS BiUT Interwa 
L=E-H 

INPUT "Hormbre d’'extrapolations: cord 


son 


":f 
"3e 


ile " 


re de la methode n = "5H 


za 


IF Hs OR CHESINTCHOS THEH FING 


GOTO 256 


® 


PRINT "Precision relative = ":; 
IHPUT P: REM ### Ph=iE-S L'4#X 


PRINT "Tableau des valeurs aPProche 
:PRIHT 

REM ### Alaorithme ##% 

DIM ICH+15 

1619 = L#CFH FCAù + FN FCB33/2 


M=z: REM ### Subhdivision de l'inter 


“alle *YX# 


3384 
344 
K 
339 
364 
30 
238 
394 
408 
414 
42 
1) 
4239 
449 


FOR K=i TO H 
H = LM: S=ÿ : I=IC1D:PRINT'K = "3 


FOR K1 = 1 TO M-1 STEF 2 

S=S + FH FCA#KI#H 5 

NEXAT K1 

ICK+15 = ICKIeEZ + S#H 

R=1 : REM ##%4 ExtraPolation ##%# 

FOR K1 = K TO 1 STEF -1 

R = 4#R 

ICK19 = ICK1+1 +0 ICK1+13-ICK1313/CRe 


PRINT ICKH12, 
NEXAT K1 
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459 PRINT: I2=IC13:IF ABSCIZ YF THEH 52 


a 

464 IF ABSCI1-129 <= ABSCF#I2: THEH PRI 
NT "integrale ICf.a.boz ";:12: EHD 

474 M=Z%M 

438 NEXT K 

499 PRINT "Precision relative "Pi" non 
obtenie. augmentez la valeur de H," 

549 PRINT "Dernieres valeurs aPProchees 
y 

519 PRINT "Ii = "511,5" [2 = ":12 

9193 PRINT"ICIL - 12394125 = ";:ABSCCII1-12 
2/12: END 

524 PRINT "I = 9 ? Changez 14 limite su 
Perieure b ": END 


DÉROULEMENT DU PROGRAMME 


Le programme calcule une valeur approchée 1(1)=—12 de l'intégrale : 
b 
I =| f(x) dx 


f étant une fonction continue, bornée, ou f donnée en forme de tableau : 


a+k(b—a) 


f=f( } k=0, …, 2N 


[a, b] désigne un intervalle borné. 
On introduit en lignes : 


100 à 110 Fonction f : DEF FN F(X)=.… 


200 à 210 A : limite inférieure; B : limite supérieure de l'inter- 
valle. 

250 Ordre N de la méthode. 

270 Précision relative P. 


77e 


Si la précision n’est pas obtenue, le programme affiche les deux 
dernières approchées 11 et 12 ainsi que l’erreur relative : 


abs (2) 


Dans ce cas-là, il faut augmenter N ou appliquer le programme 
plusieurs fois après un partage de l’intervalle (intégration par arcs). 

Si, au cours du calcul, 1(1) devient nulle, le programme s’arrête, la 
notion de précision relative n’est plus appropriée et il faut subdiviser 
l'intervalle afin d’arriver à deux intégrales différentes de zéro. 


EXEMPLES NUMÉRIQUES 


Four em Tr 

Limite inférieure 4 = 7 A 

Limite suPerieure kb = 7 1 

C4, 1 1 Intervalle 

Hombre d'extraPolations: ordre de 14 me 


thode n = TS 
Precision relative = 7 16-45 
Tableau des valeurs aFProchees: 


K = 1 
171875 
K = à 
.129159391 . 1262420483 
K = 3 
.icrereadrs . 1259424545 . 125 
Precision relative 1E-46 non chtenue, 
augmentez la valeur de N, 
Dernieres valeurs aPProchees 
I1 = .126SA0A8% 12 =  ,1275 
HOII — IZLNSI21 =  .BIGdiéééEé 


foi = KT 
Limite inférieure a = 7 Q 
Limite EUR érieure hi = + { 


CG: 1 A Intervalle 

Hombre d'extraPolations: ordre de la 
thode n = 7 4 

Precision relative = 7 1E-A5 


Tableau des valeurs aPFrochess: 


ms Hire 
fa 4j 
D 


LH 4j 
EX Le: | 
CR | 

i D 

+ 

» 

+ 

à 

Ty 

i1ù 

Et 

Fi 

ERX | 

[RE] 

[RS | 


Pi tre Di Qi re 


ef 


—} 
ci 

pes 
{ 
LE 

+ 
[ER 
cn 
Lux] 
Fo) 
Lx) 
LES 
& 
in 


125 


2 


4 
1774 . 1250443515 D ss 


0 et one ee on 





gt] 


fi = 1/S0R6.2%2inT 26 x 
Limite os a = 
Limite superieure kb = 
CA, 1.57979633 7 Inter 
Hombre d'extraPolations: ordre de 14 
thode n = 75 

Precision relative = 7 1E-095 
Tableau des valeurs afPproche 


+ 
Fe 
1. 

44 


K= 1 

2. 11980954 

K_= 

2.1516789 2.15445683 
K= 3 

2.15647145 2.156791 
2.15682845 
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me 


K = 4 

2. 156319363 
2.15651427 
K = 5 
2.12631366 
2.15651561 
2.13631362 


Precision relative 


2.13631333 
2.153631344 


2.13631366 
2.15631561 


1E-63 


augmentez 13 valeur de N. 
Dernieres valeurs aPProchees: 


I1 = 


fes 


Limit 
CG 


2.136513594 
1611 + 123-128: 


= 1/SORC, 2%sint Zac x 
Limite inferieure a = 7 

e suPerieure bh = 7 
1.57QP79655 TJ Int 
Hombre d'extraPolatior C 
thode n = 7 à 


12 = 


non obtenue, 


2. 15631562 


1.19496912E-06 


c)+co 
aa 
5 1.5 
Lersal 
ns: © 
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© 
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= ? 

Eu 15691565 13691565 
2, 13651565 2,1569156S 
24 15651565 2. 13681565 
2156921565 
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39525 
pi 1 





LR 


D'ici 


L 


EU | 







KE = 
SL ASE LELSE 
.A . HEGÉE 1254 
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. HÉT 5639 .AEGÉELATAS 
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. É je a | RÉSILIER E 
AA . RE 3P76 
Aer . 6: Arr 
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AT 
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fé = sit LAHPIE 


Limite inferieure 4 = % 1,5 
Limite superieure kb = Ÿ 


CO 1,5% 4 2 J Intervalle 
Hombre d'extrapolationus: ordre 
thode n = Ÿ 1 


Precision relative = % 1E-65 
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GÉNÉRALISATIONS 


mire 
O0 À 
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l 
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Nous reprenons la formule (7.1) et nous posons : 


n(S) v=an(s) «sin (à) 
u =sin | — }, v=sin{—}; k—sin | — 
2 2 - 2 


D'où : 


É:f* 1 
= f EE 
gl, V—-k+uV—-k-uVIi-u 


1 
(7.2) 


Nous sommes conduits à étudier le calcul d’intégrales : 





- 1 : 
dx, = , a _— 
from ee 


a 


où s présente une singularité et où la fonction f est régulière. 


= 84 = 


On peut opérer de plusieurs façons : 


1. On pose : 





(La suite (M,) est donc strictement croissante vers l'infini.) 


On calcule : 


=f" s(x) f(x) dx 


+ 


et on trouve I=I,+1,+...+.… 


2. Une ou plusieurs intégrations par parties «améliorent » les fonctions à 
intégrer. 


3. On remplace f par un polynôme approprié. 
Exemple 


Reprenons l’exemple du premier paragraphe. Soit : 
—T —T T 
Pers» a. ent et ou 


Nous trouvons pour l'écart de temps que le pendule met à passer de 
T 


à — : 


4 


‘(2 
TT/4 2 
à 
di V05 sine C) 
2 
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l'angle 


Dt= 


Tr/4 dy 
=| nn  — 06AT 


sr V2- 4 sin? (2) 
2 
Discutons encore l'intégrale généralisée : 


. 9 dx 
I= lim — 
z + — 7/2 z _— X 
2 — 4 sin” (°) 
2 


et appliquons la première méthode proposée dans ce paragraphe. Nous 
trouvons ainsi : 

, Intégrale I(f, —1,47079633, 0)= 1,40679 

, Intégrale I(f, —1,56079633, — 1,47079633)=0,30587 

, Intégrale I(f, —1,56979633, — 1,56079633) = 0,0967 

, Intégrale I(f, —1,57069633, — 1,56979633) = 0,03058 

, Intégrale I(f, —1,57078633, — 1,57069633)=9,67 E — 03 

, Intégrale I(f, —1,57079533, —1,57078633) = 3,06 E — 03 

, Intégrale I(f, —1,57079623, —1,57079533)=9,7 E — 04 


I1= 1,8537 


La seconde méthode nous fournit après une intégration de (7.2) par 
parties : 


99 Z=—1/SOR (2) 
100 DEF FN F(X)=SOR (X — Z)#(3#xXA 2 + 2x54X — 1)/((XA2—1) 
*(X+2Z))A1,5 


1=1,8541 


- 85 — 


 — 
RÉSOLUTION D'ÉQUATIONS: 
POLYNÔMES 





8. Résolution d’une équation non linéaire 
9. Polynôme caractéristique, déterminants 
10. Stabilité d’un polynôme 

11. Zéros d’un polynôme 


12. Résolution d’un système d'équations linéaires 








8. RÉSOLUTION D'UNE ÉQUATION NON LINÉAIRE 





INITIATION AU PROBLÈME 


Considérons un tronc cylindrique d’un poids spécifique de 7 N/dm* 
d’une longueur € et de diamètre d plongé dans l’eau. Jusqu’à quelle 
profondeur va-t-il disparaître sous l’eau (poids spécifique 
9,80665 N/dm*)? 


Il faut donc trouver le volume partiel V sous l’eau du tronc (voir 
ci-dessous) qui, si l’on observe la figure, s’avère être : 


V = secteur circulaire + 2 triangles 


d 2 
V=e£ (5) (7 — p<+sin @ cos p) 





La loi d’Archimède nous livre l'équation à résoudre : 
1. 
9,80665 (T7 — g+sin g cos p)—=77T <= PS sin 2@+0,89912 


d 2 
ayant éliminé le facteur commun £# (5) 5 
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L’équation équivaut (cos = x) à : 
arc cos x—xV1—x?—=0,89912, p=(1+x)d 


dont il faut calculer la solution x, 0O<x<1. 


MÉTHODE 


Il s’agit de résoudre une équation non linéaire f(x)=0. Nous 
choisissons la méthode itérative de King et Van de Vel qui améliore la 
méthode bien connue de Newton : 


f(x) 


Xn41 — Xn f'(x,) 


(x) 


f'(x) 





en cas de racines multiples où le rapport risque de prendre des 
valeurs aléatoires (1). 


Posons : 


f(x) 
fx 





U, = 


2 


et partons d’un groupe de valeurs initiales : 
Xos Mo Xi Xo — Mo Uo 
La relation itérative s’écrit donc : 


Uuy 


Mate Mn —— Xn+2— Anti Mai Un+: 
Un — Un+: 


1. R.F. King, Improving the Van de Vel root finding method. Computing 30 (1983), 
373-378. 
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VS 





L'ordre de convergence s'avère être : » lim m, la multiplicité 


de la formule en question. 


Plus l’ordre est grand, plus la convergence est rapide. 
Le programme demande la définition de la fonction f et de la dérivée 


fr 


f(x) 
30 DEF FN U(X)=(=—— 
f(x) 
Ensuite on introduit les valeurs initiales x,, m, et la précision 
envisagée. 


PROGRAMME 


14 FE XX Fesolution d'une guation nn 
nn-lingaire Fous = A KKX 

LA RENM £Æ#X Introduction de Va Fonction 
u = fsf? HAE 

Be CEF FHFC GS rer CCE fa 1 GE MEUAE 


Hi 8 EMI 
3 DEF _FN GC DEC CE 1620kH 1 dd ER +B2A TKX 
RES 


5 GEF "EN UCRieFH Fée FH Gi 

4A FEM #£#E Valeurs initiales «A, mA ### 
JA FRINT " Valeurs initiales #6, mû = " 
INPUT XG. MA: PRINT #4@:",":Ma 


6 FRIEHT " Precision p = "3: INEUTE:PARIH 
TF 

GA FENM Æ£E Alaorithme de VAH GE VEL — 
ÉING ### 


190 1 = HA — MAHFH LCHG 
119 PRINT" 1",41,M@ 

124 H=@ 

125 M1=Na 

158 MG = MA#FH LCHKBD/CFN UCKGD — FH UC 


Dr, 
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144 #a #1 

154 41 #1 

164 IF ABS: 
tion # = "3; A EHD 

174 PRINT M4, 41, MG: HE=H+1 


- MAXFH OUT 


DHIH 
mn 


£1 — #4) € P THEN PRINT" Sol 


BA IF H > 2m THEN FEIHT " Frecision no 


n chtense apres MINI jterations, 


154 GOTO 125 


EXEMPLES NUMÉRIQUES 


Tronc cylindrique 


14 FEM #EX% Resolution d'une 





gu-lineaire Fix = D ### 
EE RE EEE Totroduction de 
us pe 44 At 





Fe . GE FH LE Ge RTEC GORE 





pie LR bare SAUTE SSSR 1 — 


1 


fixisarceus x — RAGURE 1x 
Valeurs initiales x9, mA 
Frecis ion Fa = JE-BE 


(H 





Fe: we ISAFCCOS M #KkSORÉ 1"; “. 


Valeur initiale “f = 5 
rte F = 1E-HE5 
Aläorithme de HEUTOH 

.siadgeéld 
CC PACE 
RCD PAS TE 21 5 PA 
Solution x = ,.S4dZÉéÉrTÉRZ 


ONE 


F 


"+ END 


gauation mn 
la Fonction 


LR F ETHTU4 PO RAPIDES ee SÉEOReE 


—. 
; 
Qù 
LR 
= Ki 
= 
fi 


Lu 
Le 
RE 
Sn‘ ue 


h. à 
Valeurs initi “ft mA = 


Frecision F = 


1 .Sseadiiré 

E _2SCESE REA 

5 Aassréélidi 

4 ,33535634€r 

5 .s3r66es15 
Solution x = ,SSrÉ 


AR CPC Ps Cie S EU CRE PACE TT 
Maleurs initiales «A, mé = 1 
Frecision F = 1E-45 













1 .fiz9411fé 

E  . LHACSÉeSE 

5 .33r691141 

4 .333:65427 

! 3937662319 

6  .3:réd3Edd 

F dira d rez 

5. .Jasercd6 

3  .Sd6l1Azérs 

15 4645 TASSE 

11 .2133325114 

1= . 3461435966 

15 145 rars1d 

14 . 23329066 
Solution x = ,33339 

F1-3: Zero de l'ordre 4 

Maleurs initiales x, mg = -— 
Frecision F = tE— _gé 


nee 
” cbe_ 


,ABU TER 


et co x 
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cr d+ibé DA THÉ 


F 


et 1 


Fée eee ET ET dl BETETE xt 1 
Maleur initiale xf = 1 

Frecision F = 1E-64 

Alaorithme de MEUTOH 


Î -952741176 
Far sére eg 
3 "HAS TT 
4 ,r30A3E 

D mirage 

É LAarana 

Fa 44fssé 





1H 
Fes 
& 
Re 
D Lo 


. 34994381 


Cost bane Lan Las Ep Dm Le bp 
«0333204 
Lan] 


a 


++ 
£r 
8 


Lt dt CO QG) 


JT LL 'é Ne R 


[Re] 


DU Mmere £ Br ic 


LI 
à C2 Ce LOU) D Lo) Ü2 
D A9 3 LT CT EE À NA 


Ha: er BB LOMME L Où Li Ut 


EUX ne 


“Hilo ieemennk pe ke réhes ba LE 
me CS O0 €) St QC Mi Ci NU LR 
Z © ir re QU La 9 Mini 
Lt QU CN te ID QC ee CM ee M NTM 
DJS Jr Pi M St ON A 
KO CN a di O9 Ni 0 BR CRU —j Lo 


LR 'étfime ini 


A 
DEC 
ii tn 


T- 
i 


= ,.3%41 
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9. POLYNÔME CARACTÉRISTIQUE, 
DÉTERMINANT (AI — A) D'UNE MATRICE A 
ET LA MATRICE COMPLÉMENTAIRE 














INITIATION AU PROBLÈME () 


Les systèmes différentiels linéaires à coefficients constants qu’on peut 
écrire sous la forme matricielle : 


y'(x)=Ay(x)+ b(x), An, n)-matrice 
ont un grand intérêt pratique. 
On les rencontre dans la plupart des problèmes de mouvements en 
physique quand on linéarise leurs équations différentielles de la forme : 


z'=F(x z, 7) 


après avoir posé : 


s-(1 


Cette équation est d’ailleurs une conséquence du principe fondamental 
de la dynamique : 


F=ma 


Ti 


désigne la résultante des forces agissant sur la masse m; 


[SE 


désigne le vecteur accélération. 
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Prenons comme exemple le pendule pesant (voir le Chapitre 7) qui 
obéit à l’équation Z= — (£) sin z (L étant sa longueur, g constante de 


gravitation universelle). 
On introduit : 


JE 
2 Z 
On trouve : 
(2) -_ “i 

y — — sin 

2 L Yi 
et, en remplaçant sin y, par y,, on aboutit à : 
4 0 1 
Y.— 8 y 
L 


0 
L'intégration du système homogène y=A} ou : 


Yi(x)= Z ay,(x),i=1,.,n (+) 


i=1 


se ramène à un problème purement algébrique : 
en posant y = y,e*, (x) équivaut à (AI— A)y,=0. 


Ce système d'équations linéaires possède une solution non nulle si et 
seulement si À est une racine du polynôme caractéristique : 


P(A)=dét(ÀAI—A)=À"—tr AA"T'+..+(—1)" dét A 
tt A: =a,,+a,,+..+a,, 
On a donc intérêt à connaître le polynôme caractéristique. 
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METHODE 
Nous appliquons la méthode de Fadeev(*) et nous posons : 
P(A)=A"—p, À" '—p,À"T 2—...— p, 


Ensuite, nous introduisons la matrice complémentaire B(À) de AI— A, 
ce qui veut dire : 


B (À) remplit la relation : 
B(A)-(AI—A)=P(A)=dét(AI— A). 
Il en résulte : 
B(A)=B,À""'+B,À"?+...+B,_;, Bo=I 
et : 
B,=A*—p, AX1— ..— p, 1 
où : 
B,=AB,_,;—p,l, k=1,.., n—1 
Si À est inversible, on a : 


A1 Bari 
Pn 


Pa =(—1)""' dét A 
B,-_, est donc la matrice complémentaire de — A. L’algorithme de 


Fadeev nous donne le moyen de calculer simultanément les coefficients p, 


1. F.R. Gantmacher : The Theory of Matrices t. I, pp. 84-89, Chelsea Publishing 
Company, New York. 


J. Lelong-Ferrand, J.-M. Arnaudies : Cours de mathématiques t. IV, pp. 57-63, t. I, 
p- 166, Dunod, Paris, 1977. 
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et les coefficients matriciels B, de la matrice complémentaire B(À). Les 


coefficients p, sont liés aux traces des puissances A“ par les formules de 
Newton : 


tr Af=s,= X À“ k=1, 
j=1 





sh 
KP = Sx — Pi Sk—1 — + — Pk-1 Si 
On pose donc successivement : 
A, =A Pp,=tA, B,=A,—p,l 
À, = AB, p=tt—  B,=A,—p,l 
Ak 
A,= AB, ; PE B,=A,—p,1l 
An 
A,=AB,-; Ph =tr — B,=A,—p,1=0 
n 


Effectivement, on retrouve : 


AÇ=A*—p, A“ 1—...—p,_, A 


et : 


Ag = Se — Pi Ski ee — Pr Si = KPx 


ORGANIGRAMME 


P(A)=AÀ"—p; À" —p, À —..—p, 


B matrice complémentaire de — A. 


ER FR 


Entrées et Degré du polynôme n : INPUT N 
initialisation 


Coefficients a, de la matrice A par ligne : 
60, 700 


700 DATA 4,1, @, .… &in 


710 DATA @,,, 4,,, .…, a,,, etc. 


AE jÿ) <— a, B(i j) <— a; 


k=1,..,n—1 
TR=0:i=1,...,n 
TR=TR+B(i i) 
TS EL 
k 


Algorithme i=l,..,n 


200-450 B(i i)=B(à i)—P(k) 


OÙI —— résultats 


NON 


j=i1l,.,n 








B1(i)=B(i, j) (j'"e colonne) 
calcul i=1,...,n 
Axa: SS=0:€=1,...,n 
SS=SS+A(i{)B1(£€) 
B(i, j)=SS 
Résultats i=l,..,n 
500-660 Lisa 
B(i, j) - matrice complémentaire de — A 
k=1,...,n 


P(k) - coefficients p, 
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PROGRAM 


14 
15 
24 
23 
34 
39 
49 
pa) 
£a 


MONDE QD ES 


DUR EL R GDS D D St © 


C3 Fute 


TE Ep reed 
D © 


if 
Da] 


64 


ee 
# Ni 
—} 
& 


Lo cf 


UC 0) 09 Lai € 03 Où QG FOUT 


D M Di À O0) À fn in 
Ro D À © D + ER ER EG 


ME 
FEM 
REM 
REM 
REM 
REM 
REM 


FXX 
FK# 
XAX 
XX 
K#X# 
KEX 


Palsnome caracteristique#*# 
d'une matrice carree 4 *## 
Matrice complementaire ##* 
BG de -À,. HAE 
n dedre du Folynome FRA: 
À Énni-matrice HAE 


PRIMT'PEL El np ile n-15-.. pi 
FRIHT'BCL EC LEÆIT-H53 = detéL£I - Aï" 


IHPUIT 


"Ordre de la matrice n = "5 


CIM ACHHSECH HS B1CH 5, PCM 


FÜR 


1=1 


PRINT " 
FOR Ji TO H 

REAC ACT. Ji: CI, Jef TI. li 
PRIHT HCI, Jo", "; 

HET 
FRIMT 
HEXT 
FENM ##% Algorithme de FADEEY ##%# 
FOR K = 1 TO H-1 

FENM #X* Trace ### 

TR=A 
FOR 1=1 TO H 

TE = TR + EI. I: 
HET I 
PEKIETR KE 

FEM ### Calcul de 14 matrice BECK # 


FOR 


I 


TO H 
MMS Ti" Ji = "; 


1 TO MH 


ECT,I19 = BéI.It - FCK 
MEXT I 


IF K=H-1 THEH Sa@ 
FEM Æ#X Calcul de 


la matrice ACK+410 


FÜR Je TO H 
FOR I=1 TO H 
Bic lo = EI.) 
NEXT I 
FOR I=1 TO H 


SG 
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394 FOR L=1i TO H 

4n4 525255 + ACI.LoI#B16L 0) 
41 MEAT L 

426 EfI,li = S£ 

424 HEXT I 

440 HEXT I 


4SA HEXT K 

dan FENM ### Resultats #4## 

SA FRIHT 

SA PRIHT"'LA matrice complementaire EE 
CD de -ÀA ect :" 

Sc4 FOR I=1 TO H 

529 PRINT " BEM Ti" Ju = ": 
Sd FOR Jei TO H 
ESA FRINT BI. Ji" ": 
SEA HEAT 
37 so LE 
SEA MEXT 
590 RENM ### Cosfficients FCk OU EX 

684 PHIMT: PRIHT'LeS coefficients PCk, 
k=z1l,,..1n sont: 

14 FéH3=0 

éea FOR Kz1 TO H 

ESA PEN = POH5S + ACIKTHECK. 15 

6d4A PRINT " péiki"i æ MiPeKO 

ESA MEAT 

£64 END 

AAA RENM ##*# Cocefficients de 4 Par ligne 


D 


LH Li 


DÉROULEMENT DU PROGRAMME 


Partant de la matrice A, on introduit les coefficients a, à partir de la 
ligne 700 : 


700 DATA 1,2,3 première ligne 
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710 DATA 2,3,4 deuxième ligne 
720 DATA 3,45 ni" ligne 


Le programme demande le degré du polynôme n(n—=3), calcule 
ensuite les coefficients p, du polynôme caractéristique : 


P(A)= A" — p AT p A — p, 


et donne la matrice complémentaire B de — A. 
On en déduit : 


B 
ins si P, Æ0 


Pa=(—1)"! dét A 


Pour garantir la stabilité de l'algorithme, le degré n ne devrait pas 
dépasser 10. 

Des coefficients a; entiers sont préférables. 

Ensuite on peut se servir du programme Bairstow pour calculer les 


racines du polynôme. 


UNE APPLICATION 


Nous allons étudier un circuit électronique comprenant des résistances, 
des condensateurs ou des inductances, selon le schéma suivant : 


Z(s) Z(s) Z{s) 


A To. 
Vi ) Ex à) ZY (5) is Y(s) 


V5 
Sa fonction de transfert T—— présente un intérêt particulier. 
Vi 
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Nous notons les n équations du circuit : 


v (s)=(Z+Y)i(s)— Yi, (s) 


O=—Yi(s)+(Z+2Y)i,(s)— Yi,(s) 


— Yi, (s)+(Z+2Y)i,(s) 


Z Vi 
= —1 0 . 
Y Y 
2+2 1 O0 
Y 
Z 
—1  2+— —1 0 
++ 
0 0 1 2+2 0 
Y 
0 0 0 —1 O0 
—1 0 0 
2+2 —1 O0 0 
EN 
1 SRE 1 0 0 
TOY 
0 —1 2+— —1 
0 0 1 2+2 
Y 
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Y LS = n—1 
SD ee 


“ Z 
dét (15 — A) 
Y 
— 1 1. 0 
A = 2 1 
0 “1 2 


Il faut donc calculer le polynôme caractéristique de la matrice A. 
Toutes les valeurs propres sont d’ailleurs négatives et les zéros du 
polynôme également (A est définie négative). Quoiqu'il existe un 
algorithme très simple pour calculer le polynôme caractéristique d’une 
matrice tridiagonale, nous allons appliquer la méthode de Fadeev. 


Remarque 


Une matrice tridiagonale possède des éléments non nuls uniquement 
dans la diagonale principale et dans les deux diagonales adjacentes. 


Discutons encore le cas ñn =3. On trouve quand Z=R, Y =Cs: 


1 


= 
Ÿ RC?S+5R2C?s5 +6RCs+1 


et les zéros du dénominateur sont tous situés sur l’axe négatif réel. 
Admettons une tension sinusoïdale à l'entrée. On peut donc remplacer s 
par jæ et arriver à : 


1 
nn di &° + j(6RC&—R°'C'w*) 


3 


Si la partie imaginaire 6 RC © — R°C* w* disparaît, on a un déphasage 
de — Tr entre v, et v,. La pulsation correspondante w, est alors : 


. ee 
%T VRC 
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v 
et le rapport — devient : 


V2 


Pour réaliser un oscillateur, il faut donc ajouter un amplificateur qui 


La , 29 - O5 
assure un gain égal à % pour cette fréquence Ti 


c 


w 
La fréquence de coupure à —3 dB f. => 
T 





suffit à la relation : 


pe 1 


T.(jw. | 
IT: (i@.) (1+26 R°C? w2 + 13 R°C* w* + R°CS w£) 





FA 
9 
Ou : 


xX°+13x*+26x —1—0 x=RCo, 


Notre programme Baïirstow va nous livrer tous les zéros, et un calcul 
rapide nous montre déjà : 





PéLieltn-Pilten-1—... 
BCLOXCLXI-A 9 = det(L£I - A: 


Ordre de 13 matrice n = 3 


Pi 


ACL ,Js = 1 & 3 
ACe ,J) = a 3 4 
A6S ,19 = 3 4 ee) 
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La matrice comPlementaire B=B(Gi de -A est 
BC1 :Jy = —1 £ 1 
BC2 ,4ù = 2 —d ê 
BCZ ,J5 = -1 £ 1 


J 
coefficients PCKD, k=el,,,..:n sont: 
1 = 
5 = 
3 = 


EXEMPLE 2 


PÉLIELoILEEnetess sp 
BéLHILEI-A 3 = detcLET - A5 


Ordre de 14 matrice ne F5 


Hél,Js = 1 1-2 143 1-4 1-5 
ACS, Jo = 142 1458 1-4 1/5 16 
AG, Jo = 145 144 145 146 1/7 
Afd, Ji = 14 145 1-6 1/7 1-5 
ARCS = 145 146 1F 1258 1-8 


Pour cette matrice de Hilbert H., on trouve un résultat complètement 
aberrant : 


Les coefficients FCKOD, kzl,,...:n sont 
PEL = 1.78758153 

FÉES à = -,347591176 

PES 1 = 5. S35083STE-05 

péd ÿ = —1,12565954E-46 

FO 1 = 2, 1658249E-19 


valeur exacte : 3,749-10—12 
Considérons donc 2 520H, : 
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EXEMPLE 3 


FEL relié cd Fr 
EéLIHÉLEI-Hù = de : 
Or dre de la matri 
SHd 
Fi 


di = LSE , 1264 , S44 , 654 , 


T 
Pie de 


pes 
rs Fri 
D: 


dat 
Ti 
SL 
Late 
Li 
LE 
H 
Let 
RES 
EX] 
mn 
Let 
it 
ea 
‘mme * 
RES 
+ 
Fri 
ER] 


ni 
£- 
| ES 
h 
T' 
a 
CL 
Et 
RES 
pe 
fr 
Œ 
LEE] 
Ti 
it 


Éd 


Li 
Di 0n 
= [al 5 

Less 

HE 

LE 

it 

nes 

RES 

fi 

Lex | 

Let 

Ti 


14 
ci Lire 


matrice Ccomélementaire ER Hi de -H 





Se rSAdE 


PL, kel... an sont: 





Fi Ye 1 == 
FE 4 1 = 
PE d = 
pi CH. in 


et le résultat est exact. 
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EXEMPLE 4 


PéLiel ap il té nd... -Fn 
BéLYECLEIT-Hi = dets L#l - Hi 


Cid. Le la matrice nm = 6 


Ai 1 Ji = —1 Î 
AE Ji = il 2 
AL 5 = 6 Î 
Aid ls = [a 
A6 Ji = F 
AGE 1 = à É 


La. aNRise C 
Er 1 ..l1 
ECz 4.10 

ECS er 

J' 
te 0 
a î 


or 
æ 
Le 
= 


Er 4 


me fit BEN RE 


HHRHHHH 


4 
3 
2 
EL 6 Î 
Les coefficients PFLK 
pi 2 -11 


Et 


PéE = -d45 
pis 1 = Ed 
péd 5 = -Ff 
pen Où = 1 
PiE Où = —1 


EXEMPLE 5 


Ordre de 14 matrice n = 
A1 Ju = -1 1 
AéZ Ji = 1 —E 
AS 5 = 64 1 





ee fi 
mhINE 


Les coefficients FCK 
pd 5 = -5 
FÉE M = —È 
FES 1 = —-]1 


( 


. 
El 


“ù 


F 


DANS L Fr 1 


me lementaire E= EC de 


mir EEE 


com lemer taire ÉsECHi de 


mirk 
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Lex 


CE D SU LA 


! 
D 
m 


in 


es Fi QG Gi Let ns 


Qu 
1 
ne) 


UT 


las 


et 


= 
EUX 


= 
= 1 


[en 


= Fini ire i - 


EG 
à 


{ 
[AR 


he he pa 


L'exemple 4 nous livre : 


1 


A 
ST RÉCOSHIIR CS +45 RC 5 + 84R C5 +70R?C?s52+21RCsS+1 


fonction de transfert du déphaseur 6 RC; l’exemple 5 donne T;. 
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10. STABILITÉ D'UN POLYNÔME 








INITIATION AU PROBLÈME 


Considérons une boucle à verrouillage de phase d’ordre 3. Elle se 
compose de quatre parties (voir le Schéma 10.1) : 


e Un comparateur de phase CP chargé d’élaborer la différence entre les 
phases instantanées @=— œ,— @, de deux tensions v, et v,. 


Un filtre passe bas FP chargé de filtrer la composante haute fréquence 
(par exemple deux filtres actifs montés en cascade). 


Un oscillateur commandé en tension (VCO) délivrant une tension de 
sortie v,, de pulsation &, proportionnelle à sa tension d’entrée v,, sa 
pulsation centrale étant w,. 


Un intégrateur immatériel : @,= Î &, dt. 


La fonction de transfert F de notre filtre passe bas prend la forme : 


er) 
TP 


F(P=( 


Supposons la boucle au repos verrouillée sur &, et appliquons à 
l'instant 1 —0 un échelon de fréquence : 


©, = @j + AN 
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En régime de poursuite (AQ de faible valeur), le système est régi par 


l'équation différentielle linéarisée : 











oscillateurl 
contrôle 





filtre passe bas 


ViVE 
[7] 


Comparateur de phase 













v; 











Schéma 10.1. 





T; 
F=—;:K=K,K,g’(0) 
Ti 


transmittance statique de la boucle. 


La boucle ne peut reverrouiller que si toutes les solutions @ de 
l'équation tendent vers zéro quand ft tend vers l'infini. Pour cela, il faut et 


il suffit que la partie réelle des racines de son polynôme caractéristique : 


3 


1 
pA=À+PK (a+) 


soit négative. 
Passons au cas général et regardons le polynôme : 


P(x)=a,+a x+a,x+...+a,x", a #0 


Nous voulons donc savoir si tous les zéros appartiennent au demi-plan 


gauche (voir la Figure 10.2). 
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demi-plan Z 
gauche 





Figure 10.2 : Demi-plan gauche. 


MÉTHODE 


Nous appliquons l'algorithme de Routh(!) et nous posons : 


n 
bO=a,,, k=0, 1, …, nl, ni=int (5) 





ni S _ ; n—1 
b=a,,4;1, k=0, 1, …, int s 


et : 
bÜ=a,,,=0 si n est pair 


Pour i=1, .…, n, on calcule : 


ph) 

— (i+1). — Ci—1 Qi) 
de nr ATEST bf; } = d, bi D — bi) ,, 
0 


k=0, …, nl—1, 


1. P. Henrici: Applied and Computational Complex Analysis, Vol. 1, II. Wiley, New 
York, 1974-1977. 


— 111 - 


et on pose : 
CEE 0 
Le polynôme est stable si (et seulement si) toutes les fractions d, 


existent et sont positives. 
Notre exemple nous donne : 





























n =3 d, 
2FK . 
2T, 1 i=] 
T2 
2 T2 L 
TS — 2T,F?K—1 0 i=2 
” 2FK g 
T2 2 T2 | 
= TS — 0 i=3 
(2F°K) ” 2F2K 
| 








Le polynôme est donc stable si (et seulement si) : 


d,>0 ou 2F°K7,>1 


ORGANIGRAMME 


P(x)=a,+a.x+a,x +...+a,x" 


Entrées Degré du polynôme : N «— n 
40-80 (coefficients 4,, 4;, 4;, …., a, 
250 250 DATA .,., .…,.) 


k=0,1,..,N:A(k) «— a,,n 
pair : A(N+1)=0 
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N 
Initialisation N1 <—— int (;) : S «— 0 


100-135 k=0, 1,..., N1: 
B(k) «— A(2k),C(K) «— A(2k+1) 


Algorithme id; 25:23 2NS 
C(0) 
140-200 D ——, B(0)=C(0 
— y B(0=C(0) 
B (0) = 0? OUI : polynôme instable 
NON 


dE 


k=0, 1, …, N1—1: 


C(k)=D*B(k+1)—C(k+1), B(K+1)=C(k+1) 


C(N1)=0 
ET OUI ED OUI 
polynôme 
NON NON instable 
i—— i+]1 polynôme stable. 
Variables 
N Degré du polynôme. 


| N 
N1 int G) 
2 
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D Fraction décisive de l’algorithme. 

S Variable de la stabilité. 

I, K Variables de boucles locales. 
PROGRAMME 


14 FEM ##* Alacrithme de ROUTH - Stabil 
ite d'un Folsnome 44% 

24 KENM Æ## Ceure du Polunome —: coeffici 
entsz ### 

SA INPUT "'Degre n = "3H 

44 CIM ACH+15: 

54 FÜR K=A TO H 

EG FREAD ACK° 

65 MEXT 

F9 IF H = ZHINTÉN-Z5 THEN ACH4l = mA 


144 H1 = INTÉN/25% : S=<g 
114 CIM ECH15.CCN15 
124 FOR Kz4 TO Hi 

120 ECK9Y = AC2KK)N: CCKD = ACZHK+I) 

155 HEXT 

144 FOR I = 1 TO H 

154 D = CCGi/EC Gi: PRINT D: BAC à: 
IF BCA5=A THEN 594 

155 IF Dé=A THEH £=-1 

164 FÜR K=G TO H1-1 

F4 CéKi = CHBCK41T — CéK4li: EBCK+1)ecc 
K+19 

189 HEXT K 

194 CCH1=9 

2AA HEXT I 

218 IF S=-1 THEN 394 

229 PRINT "Polgynome stable, Partie reel 
le des racines < A," 

c44 FEM #X*# nti coefficients du Palunor 
e ACAD+aC lue... +aCnixtn AK 

234 DATA 1.2.4,3,5.,6 

264 EMD 
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244 FRIHT "Polgnome instable, Partie re 
elle d'une ou Flusieurs racines >=<9 " 
314 EHD 


DÉROULEMENT DU PROGRAMME 


Partant du polynôme : 
P(x)=a,+a x+..+a,x" 


on introduit les coefficients 4,, 4,, a,, …., a, en ligne 250, n étant le degré 
du polynôme. 


250 DATA .,.,.,. 
n +1 nombres 


Le programme demande le degré du polynôme : 


30 INPUT ‘Degré n=’; N 


Il calcule ensuite les fractions décisives d, de l'algorithme et les affiche sur 
l'écran. Si les n valeurs affichées sont strictement positives, le polynôme 
est stable, dans le cas contraire il est instable. 


EXEMPLES NUMÉRIQUES 


Exemple 1 


Poux = 14204 +Sutdepsts 


S 
ce 
cn 
mn 


258 DATA 
Dedre nm = 
ê 


1,5 

2.23333233 
-2,39219686 
1,95882353 
Polynonme instable, Pa 
u Plusieurs racines © 


Ca rü 


ie reelle d’une © 
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Exemple 2 


Péx) = l+tdx+6x 244 3+xt 


254 DATA 1,4,6,4,1 
Dedre n = 7 d4 
4 


5 

.8 

“e 
Polanome stable, Partie reelle des racin 
es < 


Peu) = -127+du4415%02-580 mme dt 


259 CATA -12, 4, 13: 5,3, 1 
Dédre n = 75 
-,233339333 

a 


Polynome instable, Partie reelle d'ung © 
un Plusieurs racines >=9 
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11. ZÉROS D'UN POLYNÔME 





INITIATION AU PROBLÈME 


Comme on a pu le voir précédemment, en cherchant les valeurs 
propres d’une matrice A ou la fonction de transfert d’un déphaseur nRC, 
on était ramené à résoudre l'équation algébrique : 


dét(xi—A)=p(x)=x"+a, ,x"l+..+a=0 


La résolution d’une équation différentielle linéaire d'ordre n à 
coefficients constants : 


Ha, y +... + a y + a y=0 
nous conduit d’ailleurs au même polynôme caractéristique. 
La méthode de Bairstow permettra le calcul explicite des racines d’un 
polynôme quelconque. Si le degré de p(x) devient trop grand, une 
certaine instabilité est à craindre. En ce cas-là, on attaque le problème 


des valeurs propres d’une matrice directement en se servant des méthodes 
directes de Jacobi, Householder (!). 


MÉTHODE 


Nous considérons donc une équation non linéaire : 


f(x)=0 


1. J. Legras : Méthodes et Techniques de l'analyse numérique, Dunod, Paris, 1971. 
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où f est un polynôme p de degré n : 
pix)=a,x"+a,_,x"-!+..+a,, a, £0 
Nous supposons les coefficients a, réels. 
p possède n racines réelles ou complexes et les zéros non réels sont 
conjugués deux à deux : 


p(2)=0, = %0+ jo + PR) = 0, Z0= Xo— jYo 


Pour que 2, soit une racine de multiplicité m, il faut et il suffit que : 


dd” 
p"°(t)=0 et p°"(20)Æ 0, (pw= TE) 
dx 


On peut donc décomposer p en un produit : 
p)=a, (x + pix+ qi)". (+ p,x+ q,)re (x — x, 4) (x — x) 
P du XX réels, 2(m+...+m,)+m,.;+.+m=n 
Divisons le polynôme par le trinôme x? + px + q : 
p(x)= (x + px +) p (x) +R (p, g)x+S(p, q) 
P(X)=b,X" + b, x +... + b, 
I faut choisir p et q de façon que Rp, q) et S(p, q) s’annulent 
(paramètre PS précision de solution). 
Pour trouver un couple (p, g) convenable, on applique la méthode de 
Newton. En partant des valeurs initiales p,, 4, (paramètres P, Q), on 


obtient deux nouvelles approches p,, q, et ensuite par itération p,, q, 
(paramètres NI : nombre d’itérations, PN précision) : 


CCS Cr) 0) 
q, Qn_1 R,S, —S, R, —S,R, S ÏCPn—1; An—1) 
as 
Sy = — (Pur di) etc. 
°q 
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Nous adaptons le schéma de Horner et nous trouvons : 
b,=a, 
b,-1=4,-;— pb, 
b,=a.—pb,,,—qb,,2, k=n—72,...,0 


Reste à calculer les dérivées $,, S,, R,, R,. 
En évaluant des formules de récurrence analogues : 


= b, 
Cn1= Dr-1— PC, 
Cr = by — PCxr1 — QCry2 k=n—2,...,72 
C1 — PC2— 403 
on vérifie : 
R(P, 4)= b1, Sp, 4)= b5+ pb: 
R,(p, g)= —c;, R;(p, g)= —c; 
S,(P, 4)= —c;— pc;, S,(P, 4)= —c;— pc; 
et : 


__—SR+R,S _ c;b—c;b 


R Sa — Sp RQ ci—cic; 


Ap 


S,R—R,S _c; by— cb; 


À qQ= = = 
À R,Sa—S,Ra C3—Ci1cs 


Après quelques itérations, on est ramené à une solution approchée p, q 
des équations : 


R(p,q)=0, S(p, q)=0 
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D'où : 
p(x) & pi(x) (x + px + q) 


et : 


sont deux racines du polynôme p. 
Le même procédé se déroule une autre fois. On pose : 


a,_2=b,, k=2,..,n 


n «— n—2 
et on divise (si n>2) p, par le trinôme : 


x? + px + q.. 


Finalement, on arrive à un polynôme restant de degré 2 ou 1 dont les 
zéros sont calculés directement. 

Comme toutes les racines sont obtenues d’une manière approchée, la 
précision se détériore au cours de l’algorithme (les termes R et S ne sont 
pas exactement nuls). Les dernières racines seront donc moins précises 
que les premières. , 

Une amélioration des résultats est pourtant facile. On reprend la 
méthode de Bairstow au départ en choisissant comme valeurs initiales les 


coefficients p, q, du dernier trinôme obtenus à la première phase. 


ORGANISATION DU PROGRAMME 


Initialisation 1700 REM x*x#* coefficients a(k) du polynôme 
a(n)x"+..+a(0) 
1710 DATA A(N), A(N—1), ., A(0) 
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Définition des constantes 


1030 PS—1E—4: précision de solution du système 
non linéaire R(P, Q)=0, 
S(P, Q)=0 
PN=1E—9: précision de la méthode de 
Newton 
NI= 70 : nombre d’itérations 


Introduction des variables 


1050 Degré du polynôme : N 
1130-1140 Valeurs initiales P, Q du trinôme 
X?+PX+Q 
Remarque 


Si les valeurs PS, PN et NI ne conviennent pas, l'utilisateur peut 
introduire ses propres valeurs. En cas de racines multiples, le choix de ces 
valeurs, et surtout des valeurs initiales, demande de temps en temps un 
peu de doigté. 


Résultats Affichage des racines du polynôme 
1150-1370 


Boucle du calcul 
1380-1420 


Les algorithmes 1500-1580 schéma de Horner 
1590-1690 méthode de Newton-Itération 
Résolution du système non linéaire 


R (P, Q)=S(P, Q)=0 
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PROGRAMME 


| D 


= 


( 


HE 


x | 


—, 
= 


FEN XXE Methode de EAIRSTOU Æ## 
REN XX Kacines d'un Folsnome &4%x 
RENM ŒEX 4 coefficients reels 4#* 
EE ŒÆX4# FPS: Precision de solution 
steme ÉCP,95 ef, SCp.9ù = XX 

FEM XX FH: Precision de 14 method 
HELUTOH € Tinearisationt #44 

A PE = 1E-d: PH = 16-84 2 HI=F4: HUJEHI 


pe jet pen pi 
VS LEE 1 
in 


UE ie 1 


, 
", 


1 


FE Riu more oi 


LE ee, 
1 7 
Pis tnt 


À: 


+ 
Fe 
ist 


LAS RE #4 ñ degre Qu Polunorne Fix Isa 

ner nd ic tait EX 

An REM ÆX# piece aps +9 EP lis + 

Qi + SCp,91 LEE 

AA THPUT Degree du Folynome nn = "5H 

RSA CENM Æ£LE Tutraduction des 
coefficients ##% 

A DIM AH ECHO, CH 

A PRINT MCuefficients ana n-ir. 


FH os 





a FOR K & HO TO & STEF -1 
AA FEEAD AK: PRINT "OMC, 


MENT 
RENM LXE Toteroduction des “aleurs 
iales 6,49 du trinome TEE AUX 


ee 
CA Et 


in A 


THPUT OM EE UF 

Mg ze "if 

FIHT:PRTIHT MFacines du PFolanome" 
FIHT Methode de BATESTOUM 

FO 5 ES THEM 1564 

FH 2 1 THEH PEINT -AC Gif 1 5 E 


MENT 
4 
ps 
1 
— 


ra 


 : 
ARLES 


ne = fn LT Ja Pat ee Pire 


A 

= À 

A à 

bi eu Ti 7 


121G F1 ='ACLIACEN à: O1 = AOC 
TG RE = Pise: IM = REXRE - Qc: IF IM 
SA THEH 1554 
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+ 
se Tes jen ven Hd et Eb pet jh Leb jeub hr 


UM 2e 


H° 


1310 PRINT FE + SÛR CIM 
1524 FEINT FE + SÛR & TM 





1254 TE = SRC-IM 
13ER PRINT RES" 4JENS IN 
124 PRTHT RE: jé"; In 

BEA HEN-S. IF Hem THEH END 
4 FÜR KE = 6 TU NH 

Fi = BK 





1428 COTO 1170 
4m REM EXÉ Methode de BAIRSTOU HE 
1500 HI 


151f oi . AH: Cie AH 

1524 BCH-19 = HM-12 - P#ECHO 

1538 CCH- 15 = ECH-13 - PCM 
G FOR K = H-2 


2 TO A STEF 1 
BECK se HACK PEREZ CEE HAE 
Dé ee BE <PACEKATD EC KA 
HEXT 
Céti e Céli - Eli 
REIN ÆEXX Tteration XXE 

CERN — Célia 


Lt 


il 


(EN ERR 


1 
ia! 





= 17 


LÉNEEC 1 — CCR ERA MT 

O1 HECC RS AMEEC AN — CC A MERE 1 110 

FO OCMESC EE Lun EL OPS Où AND CHESCEC A 

us 4 FPS THEN 174f 

IF ÉABSEP1I-PRÉPHU AND CAESCOA QUE 

AEH 1674 

É HI=HT4+1: TF HIGH THEN 16944 

LEA PFsPi:Gent.: GOTO 1524 

16 FE FRINT'Precision MiPSiT non obtenue 
nn isser d'autres valeurs jiuitiales" 
1684 EHD 

1690 FÉTHT"'P as de converGence aPres MH 


PP Un LT CA ER OC ET 


ROBE EU EUX 
+4 I 
Tr re 
Ho di 
4 
CNE: ES 
Fer 
+ 
E 
— 
de 
F 


QE 4 Fa ÉD dt Pal de ÉD 0 do 2 Pi CN à 
a Œ El 

Con) 

É 

"1 

< 

2 


ts 
. 
Es 
NA 


iterations, augmentez HI, ":EHG 
AR REM ### ani, afn-13,.,., 30471 ### 


lé 
1714 GCATA 1.64.4.4,0,64,4,64,-1 
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EXEMPLES NUMÉRIQUES 


Racines complexes 


Cedre du Polanome n = TT Ë 


Casffitients afntatn-1... .a64û 
1 A fi A fi 


a a G F4 


F 2 
q | 


CL 


Racines du Folahome 


Methode de EAIRSTOM 


1 
#1 
—1,.20314632E-07 +i4 1,n0444041Z7 
-1.20914632E-07 —-i* 1.414ama1z 


.PAPIBÉRFS +14 ,. 7Q718658S 
.PAPIRES7S -j# ,. 7071869885 
7407197753 +j* .74714823 
—, POPIQPPES -Jj# .PAF1RR231 


Racines multiples 


Ledre du Polunome n = 7 5 
Coefficients ani, a0n-1,.. ,a6@ù : 
224 -251 145 6 -£ 
1 
P = 7 1 
4 = 71 


Racines du Folanome 
Methode de FAIRSTOU 
.338163431 


. 27226688 
2234597401 +i* 7,7r7022464E-45 
.3245274n1 —i# 7.7F022464E-43 


-, 2452999126 


= 124 2 


Deare du Polunome n = 7 5 

Coeffi cients aénisatn-12...a( 0 : 
324 -351 1a8 € _8 
{ 

p = 7 -,.66 

4 = 7,111 


Racines du Polsncme 

Methode de BAIRSTOH 
. 3232278246 +1% 7.934S40729E-0% 
.333227846 1% 7,94544729E-G3 
.344F76A97 
. 3282312734 

-, 2980909313 


Valeurs exactes: 143, 143, 143, 143, -1/4 


1930 PS = 1E-7:PN = 1E-11: HI=99: HJ=NI 
Cedre cu Folanome fñ = Fa = 


Coefficients afniatn-19,...a0@5 
324 341 178 (ot 8 
1 
P = 7 -,66 
4 = 7 .111 


Racines du Polynome 
Methode de BAIRSTOLU 
3332026719 +1# 2.H1676662E-03 
3223326714 —i# 2,.01676662E-A5 
.326167462 
Bee 


ee ee ne et 


Le lecteur est invité à comparer le résultat avec celui obtenu par le 
programme «Résolution d’une équation non linéaire f(x) = 0 ». 

Nous reprenons l'étude des fonctions de transfert (voir programme 
Fadeev). 
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Déphaseur 6RC 


Ledre du polunome n = + 6 
Coefficients afniatn-1..,., 4045 
1 11 ds pe] F4 
21 1 
p = + 1 
q = 7 1 


Racines du Foalanome 
Methode de FAIRSTOU 
-, SALSFASATZ 
-1.29479425 

-,h55 Se 
-?2.241H27£4 
-3,.13613975 

3, PPASAT2Z 


Le régime transitoire risque d’être long à cause de la racine — 0,058 
très proche de zéro. 


Fréquence de coupure à —3 dB 


Ceare du Folynome Nn = % 6 
Coefficients afniatn-11..,.. 40 A0ù : 
1 La] 13 fi cé 
q 1 
pp = 7 1 
q = + 1 


Racines du Folunome 
Methode de BAIESTOU 


. 194086 
. 194226 
_5,2695924E-19 +i* 1.58625514 
_5.56959%4E-19 -j+ 1.58609914 
1.27591226E-19 +iX 2.24235019 
1.27591236E-19 —i* 3.24225019 
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On trouve donc : 


0,1943 
©, = 
RC 





ANNEXE 


14 FENM ZEROS D'UN FOLYHOME 4 COEFFICIEH 
TS COMPLEXES 
24 RENM ACHATHAAC HA TAC HET 4 AC, 
A6 9 3< 0 
34 REM COEFFICIENTS- REACKS, IMACK A WT 
EE 
INPUT DEGRE CU POLYHOME Hz "5H 

4 FEN Hz 

sa INPUT "PRECISION Pz ";F 

EA INPUT'HOMERE C°ITERATIOHS M1: ":ni 
74 FENM DATA REC, IMACAN,.. RE ACH 
13. IMACM-1 2. REACH, IACN: 

sf DATA l:l:-1,8,-1,-1,1.4 

199 DIM RACHD, IACHO RMCH, M2 TIME HS H 5, RC 
HO. I#CH5 

118 FOR K=@ TO H 

126 READ RACK 3, IACK 

134 HEXT 

149 RO=CRAS ET IHRAC 9 + TAC GET AE Gi aTC NE 
249 FOR Ks1 TO H 

219 RKCK I=ROKCOSE SEPIHK/H 1: LAC HK ISROESTH 
CEXPIÆK-HO 

239 HEAT 

SG REM ITERATIONS 

219 M=G 

J28 £= 

939 GOSUE 1424 

544 FOR Is1 TO H 

560 RAsRaCI O0: IAS=IACI 

S69 GOSUE ZAaQn 
579 GOSUB 54@g 
599 SeS+AESC REC TI 1 
585 RACIISRY: ICI 


… 


RY 2+ABSC TAC TI 4-13 


ï 
1 
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9394 HEXT 

GBA IF SF THEN SGG:FEM FRECISIOHN OETEH 
UE 

64 M=M+1: IF M=M1 THEH 769 

6sa GOT Sé8 

F8 FRINT"PRECISION ";P;:" MON OBTENHUE" 
719 PRINT"'ERFREUR:":S 

724 FRINT'DERNIERE AFFROCHEE :" 
739 GOTO 828 

894 REM RESULTATS 

61 PRINT"'SOLUTION:" 

829 FOR K=1i TO H 

634 PRINT "ZERO Ho "ikites "3 IMTÉRHCKO 
FP+ STEP" + LES INTEIKEK I P4 SEP 
844 NEAT 

839 EMDC 

199 REM ua DES CIFFEREMNCES 
1199 FOR J=2 TO ÿ 

111 FOR Kzl T0 1 

1128 RMC, RSR JICRHEK D 

1139 IMOJK IST XC JA-IXCK 

1159 RMCK, i=-RMC LL KO 

116Q IMCK,J5=-IMCJ KO 

1170 MEXT 

1179 RME, IS RACNO TIME, JS TACH 
1150 HEXT 

1185 RME 1, 13=RACHO: IMC 1,1 =IACH 
1139 RETURH 

can REM YALEUR CU POLYHOME 

2144 RZESRACNI:TZ=IACHS 

21169 FOR KEN TO 1 STEF -j 

L1EA RUSRA#RZ-TAXIZ+RACK-1 5 

2130 IESRAXIZ+IEARE+IACE-1 

2148 RE=RI 

2159 HEXT 

2164 FETUEH 

3RBA REM PRODUIT DES AFFROCHEES 
SAR RMSEME TI. 12 

3118 IM=IMCI. 1: 

3128 FÜR K=1 TO H-1 

3130 RMERMEÆRMOT, K+10-IME£TMET.K+10 
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3149 IM=RMEIME TI, K+1 4 IM#EMEC TI. K4+10 

3159 FM=RHM 

2160 HEXT 

3165 IF IM=Q THEH RY<RK-RZ/RM: TY=TA-TZ 
EM: RETURN 

3174 IF En=Q@ THEH RY=eRK-T2/IM:IY=T4+RZ 4 
IM:RETURH 

2175 G1=IM RM: Ge 1/6 RM+IME#C 2: = 1 CIM 
EM/01 1 

3184 RYSRK-RIZKOZ-T ZEUS 

2194 IVeL-TZHOR+RZENS 

3244 RETURH 
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12. RÉSOLUTION D'UN SYSTÈME 
D'ÉQUATIONS LINÉAIRES 








INITIATION AU PROBLÈME 


Commençons par un exemple, le Pont de Wheastone déséquilibré (voir 
la Figure 12.1). Le problème est classique en électricité. Il s’agit de 
déterminer les courants dans toutes les branches de ce circuit. La 
méthode des mailles indépendantes permet d’écrire : 


Ri+g(i—-i)+R;(i—i)=0 
g(ü—i)+R, 5 +R;(i—i;)=0 


Th+R:G—i)+Ri(i—b)=E 


D'où : 


IR +8+R:]+il[—-g]+à[—-R;]=0 
D [0]+ S(8+R, +R;]+i[—-g—R;]=0 


—R;]+Bl-R]+ [T+R;+R4J=E 


On a donc trois équations à trois inconnues. 
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Figure 12.1. 


MÉTHODE 


Nous allons résoudre le système linéaire : 
— b+Ax=0 A (n, n)— matrice, dét AZ O0 
par un procédé d’orthogonalisation. 


Pour que x soit la solution de Ax= b, il faut et il suffit que le vecteur 
(— 1, x;, .…, x,) soit orthogonal à (b,, 4,,, 4;,, …, a,), i=1,..,n. 


2 


On applique donc le procédé bien connu de Schmidt à la matrice : 


0 b, bb, … b, A (0) 
by 4, 4, … a, A(1) 
D, A A … 4, A (N) 


en commençant par la dernière ligne A (N) : 


A(N) 


(b,a, 54) Te ral 
| |A (N)| 

A(N—1) est d'abord remplacé par : 
A(N—1)—[A(N):A(N—TL)JA(N) 


A(N—1) 


—_——— : | | désignant la norme, etc. 
[A(N—1)| 


et ensuite par 


Il en résulte une première ligne A (0) orthogonale à toutes les autres. 
La solution du problème sera finalement : 


( Gi 62 2) 
X=——, ——,.., — — 


oo oo oo 


Une valeur a, très faible ou nulle indique un système incompatible où 
dét A=0. 

Le programme répète le procédé d’orthogonalisation jusqu’à ce que le 
vecteur A (0) soit suffisamment orthogonal à tous les autres (paramètre 
P : précision du procédé de Schmidt). 


Ensuite on détermine les résidus : 
r=b—Ax 
et on améliore la solution x en ajoutant celle du système : 
AXx,=r 
calculée de la même manière. 


Le procédé est réitéré (variable II) jusqu’à ce que la moyenne 
quadratique des résidus soit assez près de zéro. 
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ORGANISATION DU PROGRAMME 


Le cœur du programme est constitué par les lignes : 


1100-1250 Algorithme de l’orthonormalisation. 
1000-1060 Création de la matrice A à (N + 1) lignes et colonnes. 
1070 NI=0:NI Compteur de nombre d’orthonormalisa- 


tions. 
1100-1250 Algorithme de Schmidt 
G=N, N—1,.., 0: 
1. On calcule le produit scalaire KZ 
KZ=X A(G, ID':A(N,I) 
2. PS—KZ : on orthogonalise les vecteurs 


A (G, I) et A(R +1, I): 


A(G, I) «— A(G,I)—A(R+1,1)-PS 


pour R=N—1, …, G, et on calcule le produit 
scalaire 


KZ=X A(G,I)-A(R, I) 
3. On remplace A (G, I) par : 


A(G, I). 
SOR (KZ) 


la norme de la G“"* ligne est donc égale à 1. 


1300-1440 On contrôle si A (0, I) est suffisamment orthogonal à 
toutes les autres lignes A (G, I). 


(Les produits scalaires ne diffèrent de zéro que d’une valeur très 
faible.) 
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Remarquons finalement qu'il faut introduire l’ordre n du système : 
30 INPUT N, la précision P et la matrice A du système A x = b ainsi que 
le vecteur b en lignes 


10000 DATA... 10010 DATA... etc. 


PROGRAMME 


FEM #4 Solution d'un susteme ##4 
REM XX lineaire carre x = L ### 
FENM XX d'apres  MEMELCH. 1. ONSMAHH £ 


ER mnt mn 
Sin 


# Pl 


# 
24 INPUT "L'ordre du systeme on = MH 
4 INPUT "Precision du Procede d'orthoa 
mnalisation de SCHMICT F6 = "ir 

SA DIN ZCHHIACHME BEN IS HCH ON SMCM 
éA GOSUE £A14 

ra GOSUR SA14 

aan RENM ##* Sous-FProgrammes #4 

1AnA FEM ##E Orthonotmalisation #4 
1414 FOR K = GG TOH 

1424 FOR I = 4 TO H 

1824 AK.IN = ZCKiI2 

1449 HEAT 

1454 AK. Au = SK: ACGK M = EM 


1464 MEXT 

1874 HI=mA 

1184 FOR G = H TO A STEP -1 

1119 KZ = @: M=û 

1120 FOR R = H TO G STEP -1 

1124 PS = KZ: KZ = 4: LeRtM: Mal 

1144 FOR I = 9 TO H 

1150 &S = CGI) - ACL. IDKPS 

1169 AG: IH = SG: KZ = SOHACR. ID + KZ 
117A MEXT I 

1187 HEST F 

1240 IF KZ = @ THEN PRIHT "Det À = 97 " 
: END 

1210 PS = 1/SGRCKZ) 

1224 FOR I = 9 TO H 


FASLE 


1230 AG, I) = ACG:I EPS 
1244 HEXT I 

1254 HEXT G 

1347 PS=4 


1:14 FOR G=1 TO H 

1324 KZ=<A 

1:39 FOR I = A TO NH 

1544 KZ = ACA.IIEACG. I + KZ 

1354 MEXT 

1364 PS = ABSCKZ) + FS 

1274 HEXT 

138 KZ = PFS'H 

1594 PRINT "Orthoganalite de 14 solutio 
ni "KZ 

1448 HI= HI+1 

141 IF CKZ>PY AND CHI<4i THEN 1144 

142 IF KZ % PF THEN FRINT "Precision ": 
Fi" non obtenue ! 

1434 IF AESCACA. AN) < 1E-14 THEN PRIHT 
"Gusteme incomPatikhle 77 ":PRINT 

1448 IF ACG.4ÿ = G TMEH EMD 

1SAA S=-1/Af 4,4 

1519 REM ##X Solution *#*# 

1524 FOR K=1 TO H 

1534 ACG.K 1 = ACAKIES 

1540 HEAT 

1554 RETURN 

2AG4 REM #X# Initialisation-Introductio 
n de la matrice À ##X 

2818 FOR K=1 TO H 

2424 FOR I = 1 TO H 

PA3A4 REAC ZK. 10 

cAa4G MEXT 

2A5A MEXT 

2468 FOR K = 1 TOH 

ca?A RERC BK 

£ASA MEXT 


2434 PRINT " Systeme d'equations lineasi 
Fes : tt 

2144 FOR K = 1 TO H 

2114 PRIHT 

2124 FOR I = 1 TON 
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2134 PRINT ZK, ri" ": 
214Q HEXT 
2150 PRINT " .., ";:ECK) 
2164 MEXT 
£178 PRINT 
2134 RETURN 
SAGÛ REM ##Y Iteration -SECONC MEMBRE * 
** 
5414 II = A 
SGA FOR K = 1 TO 
SAZA SMCKS = Ecko 
SAd4a “Ki = A 
Sas HEXT 
S1An RENM ### Execution d'un Pas ##X 
5114 GOSUE 1914 
2124 PRINT " Solution aPres "illi" iter 
ations: " 
2154 FÜR KE = 1 TOH 
5144 XCKi = #CKI + ACAKO 
5154 FRINT #CK0 
F16A MEXT 
Sen4 RENM ## Calcul des resiqdus  #%% 
Sea PRINT " RBesidue: " 


H 


h 


5224 SI u 

5244 FOR K = 1 TO H 

2214 S = ECK5 

3324 FOR I = 1 TO H 

5220 5 = 5 — ZÉKITENCI) 
5244 MEXT 

SEA SMCK oi = S 

SSéA PRINT S 

SEP S1 = S1 + SES 

S3R0 HEXT 

2290 FRINT " Moucone duadratique = "is 


“H:IIF SIi-HS1E-12 THEH EHC 
444 II = II + 1 
Sd4tA IF II à 3 AND S1i/H % 1E-S THEH FAI 
HT " Susteme incomPatiblet? ":EHC 
Edca GOTO S114 
a490 REM Æ## Matrice À Par ei FLEX 
1: 


\, 
— 


14444 DATA -74,84,18,-11,-4, 
14014 CATA 14,-69,21,28.4.7 
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10424 DATA 66,-72,-5,7,1,4 

LAS DATA -12,66,-30,-27,3,-3 
14460 CATA 2.8,-7,-4,1,n 

1447A DATA di-12,4rd 01 

14440 REM #4 Second mers Li HE 
18144 CAT 1.4.4.4.7.,8 


re 


MEILLEURE APPROXIMATION EN MOYENNE 


Partons d’un support d’approximation (f,, …, t,) (ce sont les points où 
la fonction est parfaitement déterminée) et d’un ensemble de points (1,,, 


Ym) M=1, 


Nous voulons trouver un polynôme p tel que la somme : 


3 
1 M3 


CICR ETES 


1 


soit minimale. 
Pour trouver les coefficients x,, il faut résoudre le système d'équations 
linéaires : 


Ax=b 


avec : 


n 
di — D CES by = > EAN tr 


m=1 m=l 


Nous modifions donc ainsi le programme «Résolution d’un système 
linéaire. » 


24an REM #£# Initialisation-Introductio 
de la matrice À ### 
AG CIM TH VÉMI SN 


" 
2 
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QUE PRINT:PRINT "Pointe Cm, uémis 
Gi FOR M=1 TO 

15 PERD TM): YC MD: PRINTM, TM, ve M 
APD SC Met: MEXT 

FOR K=1 TO N 


GA S=4:T-0 


2050 


+ DD ID DID JDN MIAN 


T7 
Fm Fo a Po fo M FO F9 D) PO F9 


AIN RD 9 D 19 D 9990 


D 
Len 


Mmes hs pe (j 
D © Sn D 


eue O 


3 FOR M=i TON 

A SSH Mi: T=T+ CM) 

2 SEMI SCMIKTE MOD: VC M EC MATCH 
MEXAT 

Zik, LS: ROCK IST 


DIN ID ER 9 NS RIDE 


HEXT 

FOR I=2 TO H 

Gas 1 nsc 1 es { IKTC 1 

FOR Kz1 TO H-1 

ZCK: LISA K #1, I-1 1: SES K +1 
ECKALIESC KA TETE K 41 5 

MEST K:2CH,T1zS 

HET TI | 
PRINT " Systeme d'equations Tineai 


FOR K & 1 TO H 


FENM Æ£## Foints témi, umo KE&X 
69 DATA 4 


A CATA .52355987756,,.S 
326 CATA 1.570796287, 1 
59 DATA 2,617995878.,,5 

A DATA 3.14159265S4, 4 

Résultat : 
lution aPres 9  iterations 
1 = 4,7E-65 
PE ET RCE 
1 = ASE lESS 
5 = -,232211538 
5 = .N3635r08 
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MATRICE TRIDIAGONALE 


Initiation au problème 


La solution d’une équation différentielle du second ordre dépend de 
deux constantes arbitraires. Pour arriver à une solution unique on impose 
donc deux conditions supplémentaires en forme de conditions initiales ou 
de conditions aux limites, une à chaque extrémité de l'intervalle. 

Considérons le problème suivant : 


—y"+q()y=g(x), a<x<b 
Ve 
y(a)=@, y(b)=B 
y, q, g étant des fonctions à valeurs complexes. 
La méthode la plus simple pour attaquer ce problème est celle des 


différences finies (voir le Chapitre 16). 
On partage l'intervalle [a, b]en n + 1 intervalles égaux de longueur : 


(b— a) 
h = 
(n+1) 


et on pose : 





a= x Lx  <.LX=a+kh<...<x,,,=b 
x = YX), Qu = Q (XX), 8x = 8 (x). 
En remplaçant y”(x,) par la valeur approchée : 


2} + 
D°y,=2#1 pe Yk Re a 


on obtient finalement un système d’équations linéaires dont la matrice est 


tridiagonale : 
Yo — 


2 —- Yi — Ye 
conditions aux limites SR er re a = kKk=1l,..,n 


Es CO 


Yn+1 = B 
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(ou sous forme matricielle) : 


2+qh 


Remarque 


2+qh° 


0 0 


— 1 0 
2+q;h° —1 


0 Yi 
0 Y2 
0 LA 
1 Yn—1 
2+4q,h° Ya 





Une matrice tridiagonale possède des éléments non nuls uniquement 
dans la diagonale principale et les deux diagonales adjacentes. 


Méthode 


Partant de l'équation Ax= a, A tridiagonale on met la matrice A sous 


la forme : 


A=G:D 


où G est une matrice triangulaire inférieure très creuse (où la presque 
totalité des éléments est nulle) et D une matrice triangulaire supérieure : 


s, © 0 
d, 5, 0 

il, d; s; 

: .. O0 i,_, d,_; 
0. . 0 L; 


0 
0 
0 


n—1 


d* © © 0 0 
h d* 0 0 0 
e= 19 à 440 0 
0 0 0 à, d* 
1 s* 0 0 0 
0 1 s* 0 0 
D = 
Ho An LD te 
0 0 0 1 


La résolution des deux systèmes : 
Gi*=a et Dx=i* 


est immédiate. Notons alors l’algorithme : 








k=n—1,...,l: 


Remarque 
Une matrice triangulaire supérieure (inférieure) possède uniquement 


des éléments nuls en dessous (au-dessus) de la diagonale. 


Exemple numérique et programme 


Ordre de la matrice n = 
Diagonale Ke dé $, Ir dot 2 
Ciagognale Ke dé Im de : 


C9 mm 
[mr 


l'iaaonale Ke dés 5, Ir des 1 
SurdiaSonale Re sC1 3, Im sil ÿ = 2 1 
Surdiaonale Ke sCZ Im se ù = 1 —1 
Sousdiaaonale Ke 162 3,Im 162 3 = 1 —,.5 
Sousdiañonale Re 163 In 163 = GG 1 
Second membre Re 401 Im 461 0 = 3 3 
Second membre Re a6E 0, Im afé à = 4,5 3 
Second membre Re 463 3,Im 465 5 = 4 à 
Solution xe 1 1 = EH +i#i 

Solution %< Le 3 = € +i#ÿ 

Solution x 1 = 1 +ji#1 


REM #X## Sasteme d'equations #4## 
RE ### lineaires - matrice ### 

39 FENM ##*# tridiaSonale ### 

49 FENM ##* coefficients complexes ##%# 
94 RENM ### Introduction des coefficient 
z — Partie reelle, PF, imaginaire ### 

"100 INPUT "Ordre de 13 matrice n = ":M 

119 DIM RDEHD,ICDEHS, RECH=1 M, TECH=15.R 
IN, IICM5S, RACH9S, IACH) 

124 FOR K =1 TO H 

134 PRINT" Diaonale Re dé"iKi" 5, Im dc" 
Ki" = "ii: INPUT ROCKY, IDCK 5 

140 HEXT 

159 FOR K = 1 TO H-1 

164 PRINT" Surdiagonale Re sC"iKi"5, Im 
sC"iKi") = ";: INPUT RSKY ISCKOÙ 


Nr 
SNS 


=:142: 


174 MEXT 

180 FOR K=Z TO H 

194 PRIHT" Sousdiaganale Re ic"5ki"i, Im 
i6"iKk5"5 æ "ji: INPUTRICKZ IICK)D 

248 NEXT 

218 FOR Ks=i TO H 

224 PRINT" Second membre Re aC"iKi"i Ir 
a6"iKki"n æ ";: TNPUTRACK S IACKOD 

234 HEXT 

294 REM #4 Alaorithme ### 

208 S=RDE 1 DHRDE 1 24 IDE 1 Y#ID6 1 2 

310 RIG 126 RAC 1 TÉROC 1 24 I AC 1 DHIDE 1 10/85 
320 II 1226 TAC 1 SHRDE 1 RAC 1 DHIDC 1 22/8 
339 FOR K=2 TO H 

340 RSme REC Ke 1 HRDCK-14ISCK-1 KIDCK-1 5 
575 

359 1SCK=10=C ISCK-ADHRDCK-13-RSCKSTTIHID 
CK=-199S:RECK-10=RS 

360 ROCK OI=RDCKOI-RICHKHOIKRSCK-LOD4HIICKIHISS 


K=19 

379 IDKY=IDCKI-IICKOIKRSCK-A1 RICK IKISS 
K—1 

380 S=RDCKOIERODCKO+IDCKIKIDÉK 5 

399 RI=CRACKI-RICKOIKRICK-1HETICKOIKIICK- 
132275 

ci LISCIACKI-RICKIKIICK-19-IICKOIHRICK- 
19345 

419 RISK IS=RI#RODCKOI+IIKIDCKOI:IICKOSII#RD 
CKI-RIXIDCK 2 

424 NEXT 

439 FOR K=H-1 TO 1 STEP -1 

444 RICK) = RICKO-RSCKOIEKRICK4L H+ISC KIT 
ICK+109 

450 IICK) = IICKOI-RSCKIKIICK4+1D-ISCKIER 
ICK+10 

464 MEXT 

490 REM ### Result 
594 FOR K = 1 TON 
519 PRINT " Solution xe"ik;"5 = "ji RICK 
Di "+iE" IICKOD 

S29 NEXT 

#38 END 


ate ### 
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AUTRES EXEMPLES NUMÉRIQUES 


Exemple 1 


L'ordre du systeme FE 
[= es 


£ 

Precision du Procede 

de SCHMICT PF = % 1E-9 

Systeme d'equations lineaires 
—74 59 18 11 4 -8 : 

14 -63 1 28 [x] ra : 

6 7e 5 7 1 d. =: 
12 66 -58 -23 3 -3 se 

3 A 7 —à4 1 a fi 


4 -12 4 4 a 1 ns. A 


Orthogonalite de 14 solution: 
1.68920146E-45 
Orthogonalite de 14 solution: 
3.436A7987E-11 
Solution afres M  iterations: 
1.n0aa0093 
-7, 4445382E-07 
-2. 44401345 
15.9444231 
42, 09404715 
-56.0948313 
Pesiqdues : 
—1.192A929E-47 
1.782129354E-47 
5. 6746 4d42E-A48 
—-1.1920329E-A7 
—-1.49911612E-48 
1.49411617E-98 


a 


Mouenne Auadratique = 1,.473215S59E-14 


5 144" 


Exemple 2 


L'ordre du susteme n = T 6 

Précision du bracede d'orthogonalisation 
de SCHMIOT p = 7 1E-9 

Systeme d'eduations lineaires: 


-74 64 19 11 4 -6 ss 9 
14 -63 21 28 GG 7 285.0 


66 72 -5  ? 1 4 a 
12 66 -29 -23 4 -5 «… 1 
3 8 -7 -4 1 q . (a 


& .PÈE 4 4 a 1 ss. 9 


Qrthoganalite de la solution: 
“.623%4717E-08 
Orthogonalite de 13 solution: 
4,558063941E-11 

Solution apres A  iterations: 
-4û,anA11%é 

35, 444928 

155.990431 
-1934.0n259 
-3211.09892 

4996.A114 

Fesidus : 


A 
3.81469727E-ME 
7. 629394S5E-AÉ 


ss 


ouenne duadratique = 1.2126S96E-11 
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OrthoQonalite de la solution: 
4,.62F46A3E-A3 
OrthoSonalite de 14 solution: 
1.57266792E-12Z 

Solution aPres 1  iterations: 
-4n,nn9163 

35.948135 

155. A0A6E 
-174.4ñ44 
-3211.A17566 

49a6,4174Z7 

Residus: 


.B1463727E-G6 


SAN 


.42%31921E-12 
ion: 


Mouenne Auadratique = 
Orthogonalite de 14 solu 

.A613SA4SREÉ 
Qrthogonalite de 14 £& 
1.81756822E-11 
Solution apres 2  iterations: 
—-49,097169 

35, 4441536 

155.067 
-1n24.ñA444 
-3211,417852 

anse. Ai7éz 

Rezidus: 


ct FO 


lu 


Lu 


i 


Lu 
2 
2 


a 

-3.91469727E-0A 

A 

[z 

A 

a 

Mouenne Quadratique = 2.425312 1E-12 

Orthogonalite de 14 dt en 
AÉ1SEASGE 
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Orthogonalite de 14 solution: 
1.81726832E-11 
Solution apres %  iterations: 
-4f.A4n164 
25,44415 
155 1, AAAEE, 
—1A34, af44. 
—-5211.H1366 
4436.17 42 


tes i dus: 


3,81469727E-RE6 
A 
fi 
ñ 
q 


Moueniz Quadratique = 2, déSlAriE-1r 


Crthogonatite de la solution: 
HÉLSEASAEÉ 
Orthogonalite de 14 solution: 
1.8 1PSRRSTE- 11 
Solution afres À  iterations: 
44, AAA 16 
26, AAALSES 


155. GUDEET 


Mouenne duadratique = 7, 42S%1a01E-1? 
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ma1 


1 L°, 


= 
+ 
bn. 


Exemple 3 


L'ordre du syst 
Frecision cu Pr 
de SCHMICT F = 


Susteme 


1 pi 3 à 
Ds” Les 

2 c 4 : 
Le] LEA 

De] ce) " 


4 
‘4 A = 


et. 


Orthogonalite 


Qrthogonalite 
9, 11RSAA22E- 14 
Solution 4Fres 

-6,9734n143 
16, 468427 

-6,97244147 


2 ar Ad4SATIE-47 
Z,87454191E-47 
Ju LAeGaascE-Ar 
Mouenns quadra 


conditionne 
on éxkacte calcule 


de 
de 


Far le FPrograrme FADCEEY 


a 
2 


eme mn = Ÿ TS 

acede d'arthoaanalisation 

# 16-39 

tions lineaires 

LR] É 

si 4 

“ 12 
la solution:  ,.difisas 
la solution: 


15 


erations: 


te 


EX 


Exemple 4 


L'ardre du systeme n = %Ÿ Z 
Precision du Frocede d'orthogonalisation 
de SCHMIDT F = % 1E-9 

Systeme d'eduations lingaires: 

Î { Sas 


e de la solution: 


pe 17 


æ de 1a solution: 


HR 10 miû 
M Ji 

MS iù TS 
mi ju 


5 


Mosenne Auadratique = à 


UNE MÉTHODE ITÉRATIVE 


Parmi les nombreuses méthodes itératives (où l’on applique le même 
procédé plusieurs fois; voir le Chapitre 16), nous proposons celle de 
Gauss-Seidel et une généralisation que nous allons utiliser dans le 
Chapitre « Équations différentielles elliptiques ». 

Comme le programme est largement commenté, nous donnons tout de 
suite son listing. 
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Programme 


14 REM SOLUTION FAX ITERATION 
LA FEM METHOCE ADAPTÉE AUX GRANDS 
24 FENM SYSTÈMES LIHEAIFES CREUX 
44 REM ÀAX=E., À POSSECANT UHE 
SA REM DCIAGOMNALE PRINCIPALE COMIHANTE 
6A REM ACT. Ii: SOMMES ACI KM KT 
A RENM FPHFAMETRE WU CE LA METHOCE SÛF 
A FENM CSUCCESSIVE OVER COR UNDER à 
A FEM FELA#ATIOH, GHZ 
A4 INPUT "Ordre de 14 matrice À: He "; 
H 
114 IHPUT "Parametre Us "il 
124 INPUT 'Frecision du resultat Pa "OP 
154 INPUT "Hombre d'iterations HIz MiMTI 
H.I=NMI 
14 DIN ACHHO, BEMS, KCHi 
154 REM IHITIALISATIOH 
Aa FOR I1=1 TO H 
s14 FOR Kzi TO H 
cz FEHO AI.K 5 
24 HEXT:HEST 
244 FOR Iz1i TO H 
can READ FCI: #CI5i=4: FEM START 
264 HEAT 
2°a FOR Izi TO H 
284 FRIHT: FOR Kz1 TO H 
294 PRINT AI. KO. 
24 MEAT 
314 FRINT ",,." BI: 
224 HEAXT 
298 REM ALGORITHME 
ann M.I=H.J-1 
419 R=G: RENM RESIDL 
4£A FOR Izi TO H 
434 S5=f4 
44 FOR K=1 TO H 
454 SSESS+ ACIKOIKACK 
464 MEXT 
47 #=#CI 0 
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4BA ACI ISA MAC ESS-BCI DAC. I 5 

494 ReR+ABSE ARE T 50 

544 MEXT 

514 IF R<F THEN 5489 

224 IF MJ>4 THEH 469 

534 PRIHT:FRINT "Precision "5PFi"non obt. 
enue aPres MiNTi" jterations," 

339 PRINT "Solution XI) aPProchee:":C0 
TO 554 

548 FRINHT:PRINT'£Solution ACI%, Iz=1.,.,.: 
H, apres MSiMNI-Hi"iterations:" 

539 FOR I=i TO H 

269 PRINT Mkétili"i = MiIMTCHCI I /P+. SX 
F 

Sr7A NEXT 

529 ENT 

5934 REM MATRICE A PAK LIGHE 

694 CATA 4,-1,84,-1 

614 DATA -1,:4,-1.6 

ERA CATA W:-1,4,-1 

629 DATA -1,:4,-1,4 

ESA REM EBCIY. SECOHD MEMBRE 

ra CATA 1.4.4.4 
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IV 
INTERPOLATION 





13. Lissage par des fonctions de Spline 


14. Interpolation et transformée de Fourier rapide 





13. LISSAGE PAR DES FONCTIONS DE SPLINE 








INITIATION AU PROBLÈME 


Lorsqu'une fonction y est définie par un procédé expérimental ou 
numérique — on relève par exemple des valeurs à certains moments x, 
ou on résoud une équation différentielle par une méthode à pas —, on 
aimerait compléter cette fonction pour tout x par une interpolation. 

Citons un exemple et considérons un ensemble de températures 
relevées à Marseille un jour de septembre, la coordonnée x représentant 
l'heure du jour, la coordonnée y la température relevée. 


On cherche la courbe interpolante et la 
température à 14h 30 et à 21h. 





MÉTHODE 


Une certaine incertitude sur les valeurs de la fonction aux points 
d'interpolation interdit l'interpolation par des polynômes de haut degré. 
Nous employons donc des polynômes p, de troisième degré pour 
interpoler la fonction inconnue entre deux points du support d’interpola- 
tion. 
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(xs y (x)) et (CETTE y (x 4:1)) 


Px (x) = a+ be (x— x) + cu (x — x) + di (x — x) 


Le 


An 


polynômes seront définis par les relations : 
Po (Xo) = y (x) 


Pr (Xk+1) = Pr (Ak+i) = Y (ki) 








dp dp 

. Ga+)= es it EU Lin? 
d° P4 d? Pr: 

ne Ce)= ES (Xe) 


Pn—1 (x,) TT y(x,) 
Par conséquent, on arrive à une interpolation polynomiale par 
différents polynômes de bas degré en posant : 
p(x)= px), M=Ex=x,41, k=0, …, n—1 
p et ses deux dérivées existent et sont continues. 
Comme les polynômes p, ne sont pas encore déterminés d’une façon 


unique par les relations indiquées ci-dessus, on peut introduire une des 
trois conditions supplémentaires suivantes : 


1. p6(Xo) = Pr: (4,) = 0 
2. ph(x)= ph ,(x,), k=0, 1, 2 périodicité 
3. po(xo) = 4 pr-1(x,)= b 


Nous nous bornons au premier cas et nous voulons calculer les 
coefficients des polynômes p4. 


— 155 - 


Nos conditions nous mènent à un système d’équations linéaires : 


Fe Mu 20 0 &. à . 0 c g 

h, f, h, 0 0 c & 

0 h, f, h, 0 ms SU 

0 0 0 : 0 h,_; fn h,_ Cn—2 Bn—2 

0 0 0 à ? 0 h,_>2 fs Cn—1 Bn—1 
avec : 


fk = 2h + hi-;), he = Xi — Xe Yx = Y(X) 


_ 3 (Yrti — x) 3(Yx — ki) 
ES 
h, he; 


Heureusement sa matrice est triagonale, ce qui facilite la résolution 
(voir le paragraphe : « Résolution d’un système d’équations linéaires »). 
Le calcul explicite des autres coefficients ne présente aucune difficulté : 


a = Y(X)= Y%x 
An: — Yh Co = Ch = 0 


Ari (Cry +2 Cx) 


b, = 
ÿ h, 3 


FA _ CCu+1 bi) 
(3h) 
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DÉROULEMENT DU PROGRAMME 


Initialisation Il faut fournir au programme lorsqu'il le demande 
10-180 le nombre N et les points d’interpolation x(k), 
y(k); k=0,..,N, 
x(k)<x(k+1) 
Algorithme Le programme calcule les coefficients des 
200-410 polynômes : 


m=at+bh(x-x)+a(x- x) +d (xx) 


Résultats Affichage des coefficients sur l’écran 
420-570 Calcul des coordonnées d’un ou plusieurs points 
interpolés 
Partie graphique Dessin du polynôme interpolant p en haute 
580-750 résolution 
LE PROGRAMME 


1A FEM ### Lissage Par des fonctions ef 
line ##Y 

24 FRINT'Introaduction des Foints d'inte 
rpolation «<CAN CA KL ul. KE nIuCn 
" 

24 INPUT Mn = énèzgi"iH 

44 DIMAGHIECM-1 01, COH DCH-1 1 ECM-1 1, F6 
15, GÉNM-1 5 HCH-1 0, CH, VCHO 

“9 FOR K=G TO H 

ER PRINT "Point # "kit oseou = 1"; 

79 INPUT XCK3 YCKO 

28 IF K94 THEM IF XcKié= XCK-19 THEN PR 
INT'Choisissez KCkKDDKCK-19 ":CGOTO 74 

90 AKoi=Y KO): HEXT 
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149 PRINT "Toutes les valeurs jiustest 0 
SH "2 GET T$ 

114 IF T$<: "H" THEN 2984 

124 PRINT "Quel BALE voulez-vous corri 
aer7 entrez k ": CET T 

124 K=VALCTSI: PRINT "Point # "ik" mu 
148 INPUT #CK3. Ki 

15 IF K>A THEM IF XCKidz #CKe15 THEN P 
RINT'Chaisissez «Kékswke19 M:'COTO 144 
16Q IF K€H THEM IF XCKoèz KXCK+19 THEN P 
RINT'Choisissez xCkIEKCK4TS M:COTO 14 
17Q AKI=YCKOI 

18A GOTO 148% 

194 REM ### Matrice triagonale #4 

240 FOR K=A TO H-1 

21R HK) = KCK417 = KCKOI 

224 MEXT 

22A FOR K=1 TO H-1 

244 FéK5 = 2ÆCHCKO4HO Ke TL 2: GK NIET ACK+ 
12-ACK DD /HCKOD-GECACKOI-ACK 1 D IHCK ST 0 

£SA MEXT 

264 REM ### Resolution du systeme linex 
ire #X*#*# 

278 C1 = Gé12/FC 1 

22A FOR Kz2 TO H-1 

294 ECK- D = HOK=13/FCK-10 

20@ FK3 = FCK9Y — HOK-1Y#ECK-10 

319 CéKD =CGCK Si — HCKe1 TÉCCK-A1 INF KO 
22A HEXT 

324 FOR K = H-2 TO 1 STEP -1 

24Q CeKi = CEKD — ECKOIKCCK+10 

259 MEXT 

36A RENM ### Coefficients des Palanome 
#XY 

37A CCG = A: CON = 

2e9 FOR K=Q TO H-1 

394 ECKS = CACK4LOI-ACKOID/HCKOD CCOCK41 
+ 2HCCKOMEHCK T3 

488 DK = CCCK41N — COK INC 3KMCK 
41A HEXT 


OI 
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A2 REM #£# Affichañe des coefficients 
##+ 

43A PRIHT"'Polanomes PKCxXI, k=G,.., , n-1 
44 PRINT"'ak+hkém-xk D4cké sms 024 dk use 
Kia " 

454 FOR K = @G TO N-1 

464 PRINT "k = "Ki" sk = MiMCKOUI" x 
k+li em Mie K41 0 

47Q PRINT AK2BCKICKI DK) 

484 HEXT 

4934 FRIHT:PRINT M Calcul des coordonn 
es d'un Point IinterPole " 

AA IMPUT “Valeur de 14 variable x = ": 
Le 


S1Q IF H<KCAS OR K+ KCNS THEN PRINT "# 
en dehors de l'intervalle ":COTO Sam 
sea K=g 

SA IF ArkCK4+19 THEN KzK+1: GOTO 534 
25 Az AC 

4 PRIMNT'PEM HACK 5" = MS CCDCKIKA4CC 
KA IKHABRCKOIUIEX +ACKO 

SSA PRIHT'"'D'autres Pointe 7 O/M ": GET 
T$ 

SA IF TH<e "M" THEH SGA 

va REM ### Courke obtenue Par interphol 
ation *#*# 

529 PRIHT "Dessin sur l'ecran 7 Q/H ": 
GET T$: IF T$="N" THEN ENC 

594 INPUT "'Ordonmeesz extremes 4 min, 4 
mas ":%1.%2 

6AA IF 1 ?= Ÿ2 THEN PRINT "O4 min € 4 
max ! ": GOTO 594 

ERA HIRES:REM Ecran Haute resolution 
é24 K = A 
ESA FOR I = A TO 233 
EAD x = HAN + CHCMI-HCAIIKI/239 

CSD IF HèkCK4+19 AHD K<H-1 THEH K=k+1: 
GOTO 65 

Sn = 4H 


PRES | 


EE = CÉDÉKIKE + COKOIIKK + BCKODEX + 
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E7A = INT CONZ-NU4191/002 104,50 
ESA IF Y<4 THEN Y=f 

ESA IF 2191 THEN =191 

FAaRn CURSET I.%,1 

71R MEAT 

ren CURSET 4.4.1 

724 CRAU 9,191.1 

74 CEA 239.08,1 

SA POKESE.S6:REldPermet d'ecrire dans 1 
és 2 lignes en bas de l'ecran Haute Res 
774 PRINT'Enfoncer une touche Pour cont 
inner": 

724 GET A$ 

244 TEXT 

214 EHD 


QUELQUES EXEMPLES NUMÉRIQUES 


Courbe de températures 


Introduction des Points d'interPolatia 
n kB u A) LIN)... .kénIsaeni 
ñh = énds2iT 8 

Poivit # @ %x,4 = 7 Q6,16.5 

Point # 1 x.1 = 7 48.18.39 

Point # 2 %;u = 7 

Point # 23 x:u = 7? 12.28,4 

Point # 4 %x,4 = 7 14,29.6 

Point # 5 %.1 # 7 16,28 

Point # 6 x,4 = 7 18,22,1 

Point # 7 x,u = 7 24,19.4 

Point # 58 x.1 = 7 22,17.3 


‘Toutes les valeurs justes? O/M 

Quel Point voulez-vous corriger? entrez 
k 

Point # 5 %x,1 = 7 

? 16,28.7 

Toutes les valeurs justes? D/M 
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Polunomesz PK, keD, on 1 : 
ak+bké ssl eh ssh 1024dhé sue 373 
k= A ski = 6 x(k+15 = 9 


16,5 . 2282196774 a 
.A787253063 

k = 1 KkKI = 8 xCk+19 = 14 
18.9 1.22978645 . 4723239838 

8811165315 

k=e Z Kiki = 19 xCk+19 = 12 
23.8 2. 14374742 -. 2143593315 
-, 19425718 

k = 3 x(k) = 12 KCk#+1 3 = 14 
28.4 1.43722386 -. 63994243 
. 11864525 

k = 4 KÉkKD = 14 xCk+H1) = 16 
29,6 . 207337 144 .-n239694719 
—, 1732823822 

k = SO xékKie 16 xCk+1) = 18 
28.7  -1.80863243 -1.9349738€ 
. 2676539837 

kæ 6 xék3 = 18 xCk+19 = 20 
23.1 2. 12474743 . 74926362 
—,. 6977763254 

kæ 7 KxCkÿ = 20 xCk412 = 22 
19.4  -1,25835788 . 15626641 
—. 26044735 


Calcul des courdonnees d’un Point inter 
Fole 

Valeur de 14 variable x = 7 
7 14,5 

PE 14,5% 9 =  29.6366929 
D'autres Points 7 DH 
Valeur de 14 variable x = 7 
7 21 

PC 21 ) = 18.271865 
D'autres Pointe 7 ON 
Déssin sur l'ecran ? O/H 


un 
ic 
v] 


Ordonnees extremes 4 min. 9 max: 16.5. 
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°C 30: 








+ 
10 


Introduction des P 
n xlAÿ.g ai wkel1),: 
he én2e397 3 


in 
La 


cints d’ 
CLS D PRES A 


[4 


Paint # 
Point # 
Point # 
Paint # 
Toutes les 
Polunomes F 


D forms 


Kad 
Ha 
“a 
Ka 
valeurs 
k£wæls K 


HE 


# 


nl 
4 
! 


2 0 + 
UN TU 
+R 


jJué 


eG,, 


ak+bkCw=xk >+cke est 724 


mn 
= 


Fe A ee 
{, AGEEGEE? 


1 ue 


\ 4 


q 


«C++ 


3, ACEGEEET 


3 1 


#ek4IÙ = 
. 1323333233 


#6 Kk4 
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eZ 
dc 
12 


is 


12 


“4 
” 
€ 
a 


Len] 


3 


interpolstic 
nie 


Q“M 
1 


uk 32 
= 1 


e 
as 


3 


_ 


8 


Calcul des coordonneesz 


Pale 

Valeur de 1a variable x 
?v 2 

pe 3 = à 

D'autres Points 7 0/H 
Valeur de 14 variable x 
1.0 

PC, 7 =. 1,4 

D'autres Points 7 ON 
Valeur de 1a variable % 
7 2,8 

PE 2,8 3 = 3,16€ 
D'autres Points 7 OH 
Dessin sur l'ecran ? /H 
Crdonness extreme 
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=. 


d'un Point inter 


& 9 min. 4 max: 9, 4 





14. INTERPOLATION ET TRANSFORMÉE DE FOURIER 
© RAPIDE (FFT) 








INITIATION AU PROBLÈME 


Les séries de Fourier ainsi que la transformée de Fourier sont depuis 
longtemps des outils très puissants pour traiter les fonctions périodiques 
et non périodiques qui décrivent les phénomènes de la physique et de 
l'électronique. Citons l'amplitude d’un son en fonction du temps ou plus 
généralement le traitement de signaux et le calcul des spectres. 

Le théorème de Fourier dit qu’à toute fonction g intégrable et 
périodique peut se substituer une série de fonctions sinusoïdales dont les 
pulsations sont multiples de la pulsation fondamentale : 


27 
&=27Tf=— 
f T 
AG ee : 
g()E—+ Ÿ (A, cos nwt+B, sin nwt) 


n 


Les coefficients sont des intégrales définies : 


2 2-7 
Ar] g(t) cos nwtdt, B,=K | g(t) sin nwtdt 


0 


A 
et un désigne la valeur moyenne de la fonction g. Les fréquences nf et les 


coefficients À, B, caractérisent donc le spectre de la fonction g (voir la 
Figure 14.1). 


L 64: 


Dans le cas de fonctions non périodiques, la série de Fourier est 
remplacée par la transformée de même nom : 


spæ/] g(t)e-2rif dt 
et le spectre perd son caractère discontinu. 


On peut reconstituer la fonction g de son spectre S(f) par la 
transformée inverse : 


eo=] S(f)e Ti" df 


(si g est continue et si les intégrales : 


[__lela [ Islr 


existent). 


Comme les intégrales sont généralement difficiles à évaluer et que 
souvent on ne connaît de la fonction g que quelques valeurs, on les 
remplace par des sommes en nombre fini de la forme : 


N—1 
S(k)= ZX gee TEEN, k=0, 1,..., N—1 


{=0 


MÉTHODE 


Soit f une fonction périodique de période 2 7. 
Prenons une base d’interpolation de fonctions trigonométriques : 


1, cos x, sin x, cos 2 x, Sin 2 x, .…., COS nX, Sin ñX, … 
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Nous voulons interpoler la fonction f, dont on ne connaît les valeurs 
que dans les points x,, par l’expression trigonométrique : 


M 
È A, cos kx+B, sin kx, N=2M+1 


k=1 


A5 
AS og 5 


A6, MT! ; Am 
p(x)=— + Y A, cos kx+B, sin Re eos Mr N=2M 


k=1 
et nous admettons que les points d’interpolation : 


(x f(x)), k=0, ES | N—1 


sont équipartis sur [0,2 7] : 





x = » k=0,.., N—1 
On peut considérablement simplifier le problème en posant : 


ek% = cos kx + j sin kx (formule de Moivre) 


et: 


p(xX)= bi+ be" + be * +... + be, ex 


p polynôme de phase, b, coefficient complexe. 
On vérifie facilement : 


A A4 — jB, A4 B 
= b, = k - bi pen ES _- x) 


bu A = bi + x x Bi = jh — by 4), k=1,.… 


tp) 


A6 = 


Les polynômes de phase : 


p,(x)= >» bel, PR-1—P 
k 


=0 
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ont les qualités suivantes : 


a. p(x)=f.=f(x), k=0,1,.., N—1, f, complexe 


si et seulement si : 


1 N—1 | 
b=— À fre-2TikN k=0, 1, N—1 (14.1) 
N #=0 


b. (14.1) a pour conséquence que p, est la meilleure approximation de f 
en moyenne quadratique. 
p, est déterminé d’une façon unique par cette propriété. 


Nous nous limitons maintenant au cas N = 2". Cooley et Tukey (1965) 
ont généralisé une méthode de Gauss pour évaluer rapidement les 
expressions (14.1) avec O(N log N) opérations au lieu de O(N?), 
nécessité absolue pour effectuer le calcul de la transformée en temps réel. 


La transformée discrète est décomposée successivement en deux, 
, ; NeNn ne, s : | 
quatre, huit, etc. transformées sur Sr points jusqu’à deux points 


(voir par exemple Microsystèmes septembre-octobre 1981, pp. 155- 
159). 


Ne pouvant citer tous les détails, nous nous contenterons de dire qu'il 
faut introduire l’application T qui inverse la représentation binaire d’un 
nombre entier k : 


k=a,+a2+a,2?+...+a, 2" , a =0 où 1 
T(k)=a,_;+a,_,2+...+a,2" "7! 
on initialise donc par : 


b(rT(k))=f, k=0,.., 2"—1 
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exécution, 
coefficients du polynôme de phase multipliés par N 


Enr me se ES ss € Sn Se 
D. Re EN LE 


a F0: 2er 1 
r=0, 2", 2-2, (<N—1) 
u:=b(r+i), v:=b(r+i+2"-I).E 
b(r+i)j:=u+v, b(r+i+2" th =u—v 


—2T)j 
E <— E exp ( pr ) 





Après on trouve dans le tableau b(k}), k—0, …, N—1, les 


— mi 
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ORGANIGRAMME 


Entrée des données numériques : 


60 N nombre de points d’interpolation 
N=2\ 
110 RB(K) partie réelle de f 


IB(K) partie imaginaire de f : 





Kx*x2*PI] 
RB(K)+ jIB(K)=/( =) 


Inversion binaire : 


130 à 185 T(K)=L, T(L)=K 


Algorithme de Cooley-Tukey : 
200 à 320 Calcul des coefficients b(k) 


Affichage des résultats : 
400 à 420 b(k) et |b(k)f 


Dessin du spectre : 


430 à 470 


PROGRAMME 


19 RE ### Transformation de FOURIER -F 
FT #X# 
PA PRINT'hC x ÿsbDtblexhé ixi+bzexP(2ix ie, 
+ HEH-13exkPCEN-1 ik 
39 PRINT'&xCk)S æ ZPI KM, CxCkU, FÉKEk 332 
Points d'interpolation, k=@,1,,.,.,N-1" 
48 PRINT'Me2%n, calcul des coefficients 
hék 9" 
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59 FEM ##*X Introduction des valeurs ( 
complexes) FOxCk IIIe OKI KXKXK 

6A INPUT "Hombre de Points d'interpolat 
{on Nz2%n = M"; 

78 H1 = INTS LOGCNS-LOGE 234,55: DIMACNI 


) 

84 DIM RBCH-13, IBCN-15 

94 MI=N/2:L=A 

199 FÜR K=Q TO N-1 

195 xX=K#2/H 

119 PRINT "Re FC" 4; UEPI TN , Ir fé" ai" 
#PIS = ": 

129 INPUT RRCK 3, IBCK 0 

134 IF K<L THEN 174 

149 RB=RBEL 3: IB=IBCL) 

154 REBCLOISRECK 1: IBLI=IBCK 

164 RBCKOI=RB: IBCKISIE 

174 .J=M1 

17% IF Jé=L THEH Lel-J:JsJ/2:GO0OTQ0 175 
184 Lel+J 

183 MEXT 

Po ### Alsorithme de COÛLEY - TUKE 
cu 4 

204 Mi=1: FOR M=i TO M1 

214 A=COSC PI M1 2 : B=-SINCFI/M1 02 FR 
E=1: I1E=Q 

224 FÜR J=<A TO M1-1 

234 FOR R = @ TO N-1 STEF 2#11 

249 RA=R+.I:RM=RY+MI 

254 RU = RBCRQS: IL = IECRGS 

264 RY RECRMI£RE - IBCRMOI#IE 

27a I" RECRMIXIE +TECRMS#FE 

284 RBCROY = RIHEW: IBCRGD = IL + IV 
294 RBCRMS = RU-RY :IBCRMI = IU - I 
388 HEXT 

395 RS=AYRE-B#IE : IE=H#IE+EXRE : RESRS 
214 HEXT 

315 MI=2*%M1 

324 HEXT 

399 PRINT "Coefficients HCkO: " 

484 FÜR K = @ TO H-1 


n tt 
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414 PRIHT " KHEMS Kit = Mi RBCK MH" IE 
"5 IBCKM-H; 

415 IBSKOI=RECKOIERECE DC HÆHO+TBCKOUKIECK à 
SÉMEHi: PRINT" HE KINUISE = M: TBCKS 
42A HEXT K 

434 PRINT " CRAPHE Cl SPECTRE CE FLIISSA 
HCE CU SICHAL " 

449 FO KzA TO H-1 

454 PRINHTK : SPCCSHINTEMESORE TECH 4,510: 
LE A0 

46 MEXT 

47A EHD 


DÉROULEMENT DU PROGRAMME 


Le programme calcule les coefficients du polynôme de phase b(k). On 
introduit en ligne 60 le nombre N de points d’interpolation, N = 2", en 
| h ; …. .. .{2Tk 2 TTk 

ligne 120 les valeurs réelles et imaginaires Ré (+) Im (+) 


On place alternativement ces valeurs dans un tableau où l’on se sert 
des instructions READ, DATA. 

Les coefficients b, sont liés aux coefficients réels A, et B, du 
polynôme trigonométrique d’interpolation par les relations : 


A0=2 bo, Ax = D + bx_x Br = j(bx — bx_4) 


Pour calculer la transformée inverse, on emploie le même programme 
en modifiant : 


110 PRINT Re b(k), Im b(k)”; 


PI PI 
210 A=COS (x) : B=SIN (Go) 
M1 M1 


et en supprimant la division par N en ligne 410 : 
410 PRINT ‘’f(2pi”; K; ‘/”; N;°)="”; RB(K); ‘’+jx*x"”; IB(K) 
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EXEMPLES NUMÉRIQUES 


Exemple 1 


FéskishA+hlexkpé ik itheexPé Ze. + Hé Hi 5 
eaxPCCH- 1 x 

ski = ZPI KM, CkCk M FER K 18 Points d 
“interPalatio.s 

k=ñ,1,,..,.14N-1, Hegfn, calcul des coeffic 
ients bék 





FO ÆPFIS , Im FC #FIs = Es] 
fie API) » Im fC.25 XPI2 = , 45 
Pose API ss ne PÉCS-AFT: = 45 
Fer #PID à Im F6. FS 4AFI = F8 
FA ÆFT Im SCA FT: = 1 

FOIS ÆPI ST . Im Fél.ES #FI5 = 

FEA.S PIS, Im FC1.S #FIS = 1.5 
fol, > Im FC1.75 #PID = 


Pal 
+ 
Tr 


_ 


SE 


Ti 


DE CONS CRE CR CSL La im NE 


UM Jin 
à 1 


FT ra 
F 


Æ 
efficients s 
ét = AFS “ 7 Hé NI m TEE 
Ed nes 125 + ESA cl it & ,Igf 
DéE nee lES +4 JE 18S ikbée pit æ Gi 
HS De —,185 + JE Se his it © Gi 
Ed Me AL + JE ht d item ,AlÉ 
HES se AZ + IE OS ls it = is 
ie ne ILES + JE TES bé HUE & , dl 
HET ne 18 + JA GE pr tite = ,iOr 


GÉAFHE QU SPECTRE DE PUISSANCE CU STGHAL 
fi # 
# 


Pire 


HT LA He 
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Hombre de Points d'interpolation Heëtn & & 


Pre es Lt es L'Oeil de + Ds EG 1 





Re FE APT; Im FE ÆFIS = à Fi 

Re PO ES #PIN : Im FÉES #PIS = 2, GES 
re Peso APIs Tr PSS API) = E 

Re Fée ÆPIs , Ir F6, 75 #Plo = 3.525 
Ke FL #PI A, Im fé FI = Ë Us 
Fe Files #PI A . In FOl.ES #Fli = 
Re FEI. S ÆFT A , Im FOl,S #Fls = € 
Fe fFé1.75 ÆPII , Im F1. Fi HP = -2,SéS 





Coefficients Hé 1: 

LD = G + 48 RCE site = G 

El ne M + IE, il it = , ES 
He ne A + jE ES te it & GE 
Hé 5 = M + iE-A,S ICE nie æ PES 
En 4 OÙ oO + EG  ibéd TE & 6 

HÉS = A + J41.5 FE à HUE Er 
Hé ne Et + JE ES bé é DIE & QE 
Hé Us M + JE be r NE = LES 


GEAPHE CII SPECTRE GE PUTESANCE CL STGHAL 
# 
4 


aù Pres it 


HDI EL 
# 
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Li Li: 


Exemple 3 


PE ik tb dekPe Let + Ne 12 


Le 


Um nn x Le 


sf 


EM 


Le 
D 
k#H, CKkek fé mCk 358 Points d 
i 


CO D D EL EX 
he Te 


1: He en. calcul des coeffic 


me 


Hombre de Foints d'interFolation Heës & 16 
Fées nl es Los D is Lol 4x4, ei ré BK 


Fe fi ÆPIT: ., In fé@ #FIi = 

Re FO AE #PID . Im fé.1ES FI 

Fe FÉES ÆFITS , Im FÉES PIS = 

Fe POSTES #FIS : Im FES SFE FI 

Fe Fi,S ÆPIS , Im FES PI = -£ 










Erin 


LR Pt Gi re 
Fa 


re FOLEES #FIS , Im fé.6ES #Fli = -2,578 4 
Fe fC.75 PIS : Im f6.7S #PlD = 3.328 : 
(EF = 4,51 [s) 


LABS PIN Im FRS ÆFI 
FEI FI . Im FE #FIi = 
ne FOLIES ÆFIS , Im FÉL,122 #FI: 
Fe FOOT. ES ÆPFI M, Im FÉL,ES #FIi = 
Fe f1. 375 ÆFT OU , Im FÉA.SFS #FI 5 
Fe fe 1.5 £PI SH Im FÉ1.S FI = € 
5 #PFITS : Im ee #PI 2, 
#PIS , Im FÉ1 #Pli = -2,53È 
#PI . Im Fe. &7S HP = -1 
HE 
i+A 
It, 
it 
Ji, 5 
J#1 
i#g 
ik, ES 
RE 
#6 
itA 


3 
Lux] 
rer 


HinBii 
NE 


{ 
RCE 


in Fa 


NN Er Q 
EU 


DO 


LD 
En 


fl 
Fi 


C1. 
FT 
le 


Ë 
F 
& 
1e 


mi 
+ 4 + # 
pas ST 
RAGE 
in 


4 
Le 
Ds JU Et 


nu m 
mn mie he be LD Ou M Un La ne © 
ER Or 1 it 


EUHENMEUNMN E 


HHE+++++ ++ 


_ 
ul 


HUF GR Be RO NIP a rire ete 


fé CR n = A + i#.: 

Hé di = EH + es ati : 
HÉAE 9 = 4 + j$-1 Hi 
HAS 1 = À + #1, 1615 
Ed 1 = 6 + j£6 bclid 
Hé 1 = 4 + ES is 1° 
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Es 
Lt 


GRAPHE CL SPECTRE GE PUISSAHNCE OÙ SIGNAL 


# 


DUJDNDNRUMEE 


me pe pe Ur 
2 FO 


En 
Bio 


Le 
on 








Spectre 


Figures 14-1 
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sont HON-1exPé CNT nix 

Points d'interPolation. 
£ des coefficients HCkû 
Hombre de Points d'interpolation Neï%n = 178 






max Cki = 1E76<H 


f(x) = x/7 Spectre | b (k) | 














Fixo = sine 18#xi + bruit Féxi = sin 10. SÉx) + bruit 


max 1EkIi = 63.9264456/N max RCE = 43,33546955 /N 


Spectres | b(k)| : 


0,499 


0 











Figures 14.2 
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V 
ÉQUATIONS 
DIFFÉRENTIELLES 





15. Équations différentielles du second ordre non 
linéaires 


16. Équations différentielles aux dérivées partielles 





15. ÉQUATIONS DIFFÉRENTIELLES 
DU SECOND ORDRE NON LINÉAIRES 








INITIATION AU PROBLÈME 


Il existe en électronique une classe de montages appelés oscillateurs et 
qui jouent un rôle important. Ces montages ont pour but de fournir sur 
une de leurs bornes un signal sinusoïdal dont on maîtrise en général la 
fréquence. 

I s’agit de systèmes bouclés, dont la réaction est positive, et qui 
obéissent au critère de Barkhausen : le gain de boucle doit être égal à 
+ 1 pour une certaine fréquence, celle de l’oscillation justement. Or, il 
est difficile de réaliser une telle condition de façon parfaite. On s'arrange 
d’ailleurs en général pour que ce gain de boucle soit légèrement supérieur 
à +1 et que l’on soit sûr, par conséquent, qu’il y a oscillation. 

Par ailleurs, l’amplificateur (nécessaire) n’a pas une fonction de 
transfert parfaite du type v, —Kv,, mais il présente des paliers de 
saturation : + Vsat. 


Ampli- 
ficateur 








Figure 15.1. 
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Figure 15.2 : Fonction de transfert d'un amplificateur réel. 


La partie linéaire de cette fonction correspond au gain de l’amplifica- 
teur. Il est possible de rendre compte de l'allure de cette caractéristique 


Va ) 
Vo 


ce qui est d’ailleurs proche de la réalité, le dessin ci-dessus n'étant en fait 
qu’une modélisation de l’amplificateur. 

Cette relation, jointe à l’équation différentielle du Pont de Wien (voir 
encadré), fournit l'équation régissant le fonctionnement d’un oscillateur à 
Pont de Wien. 

L’équation devient alors (v, =: v) : 





Vi = Via tanh ( 


dv 3dv Vs 1 dv 1 
sn —+— v=0 


dt Tdt TV» En 
cosh 





20 
et on a une équation différentielle de Liénard si : 


V 
34e e>0 


20 





Cette équation possède une solution périodique stable v,, (sans prendre 
en considération les translations de temps t —+> 1+ Dt) dont toutes les 
autres solutions s’approchent (voir la Figure 15.3). 


Nous allons donc étudier le traitement numérique de cette équation. 
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Équation différentielle du Pont de Wien 


R "ee 


i(t) Es | 


vi (t) 


v; =0 


d°? v, a ( 1 1 ) 1 dv, 


—— a ——— + — és, ER 
dé a \RO'RC RC) RC & ‘RRCC 


On pose en général, pour simplifier : 








É> 


Courbe UC x à —- transe 


Figure 15.3. 


Évoquons encore un autre aspect : 
Une des équations fondamentales de la physique est la célèbre relation 
de Newton : 


F= ma 


(on masse d’un corps, F force qui peut résulter d’un potentiel lié à la 





dy 
masse; a accélération, a = Pre ) 
X 


y(x) donne la position du corps, y(x,) sa position à l'instant x,, y’(x,) sa 
vitesse initiale, x le temps. 


Comme F est une fonction des variables x, y, y’, notre problème 
s'exprime au moyen d’une équation différentielle du second ordre : 


Nr — 


F À à 
'=—=f(x y, y) 
m 
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et deux conditions initiales : 
Y (x) = yo Y'(xo) = 4 

ou si l’on pose : 
Y=}Y, ÿ'=}, par un système 
1G)= y»), Y1 (0) = Yor 


(x)= f(x (x), »(x)), y: (%0) = Yo 


MÉTHODE 


Tracé d’une solution de l'équation y"—= f(x, y, y”). 


Nous discutons d’abord le problème à valeurs initiales du second 
ordre : 


Y'= f(x y, y), y(a)= y, y'(a)=w 


et nous voulons tracer sur l’écran la courbe solution y(x) ou sa dérivée 

y'(x), aZ= x = b, et afficher les valeurs approchées dans un tableau. 
Dans ce but, nous partageons le segment [a, b] en 200 intervalles 

partiels égaux et nous utilisons sur chacun d’eux une formule de 

Runge-Kutta d'ordre 5 (erreur locale = O(n°)). 

(b— a) 


Partons donc du point initial x,, yo, w, et posons h — 700 





L1—hf(x, Yo, Uo) 


h hu, L1 
Fr hrs ES pee) y mis 
L2=hf(x, y, y’) 


hu. Lih L2 


À 


h 
nr PES Er y ut 


1484 


L3=hf(x, y, y’) 
E2h: : 
xX=X+h, Y= Yo +hu TE y'= us + L3 
L4= hf(x, y, y) 


et finalement le procédé repart à partir de : 


h(L1+L2+L3) 


Xo — Xo+h, Yo <— Yo+ hu, + 6 


(L1+2L2+2L3+LA4) 
6 


UGo — Uo 


Si l'intervalle [a, b] est important, un partage est recommandé. 


Système d'équations différentielles d'ordre 1 d'après 
Fehlberg (*) 


Nous allons étudier le système d’équations différentielles : 


IA 


vi fe (à as Vas es Mn) k=1,2,.., n, aZ=x<=b 


ou en forme vectorielle : 


u—=y =f(x, y) 


Soit x, y, le point initial (x,, 1) de la courbe solution. L'équation 
différentielle nous donne alors la pente f(x,, y,) de la tangente à la 
solution. 


* E. Fehlberg : Klassische Runge-Kutta Formeln vierter und niedrigerer Ordnung mit 
Schrittweitenkontrolle und ihre Anwendung auf Wärmeleitungsprobleme. Computing 6 
(1970), pp. 61-71. 


J. Legras : Méthodes et Techniques de l'analyse numérique, Dunod, Paris, 1971. 
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L’approximation consiste à suivre la tangente au lieu de la courbe sur 
un petit intervalle [x,, x, + h], h pas d'intégration. 
Nous posons : 


L= hf(x,, yo) 


Pour améliorer l’approximation, nous reculons et nous recommençons 
le même procédé à partir du point : 


X+ah, pm +b;L 
Cela nous donne : 
M=hf(x,+a; h, y, + b;,L) 
et le prochain point de départ sera : 
X+a@h, y+b L+b;M 


etc. La valeur approchée y,+Dy de la courbe solution devient 
finalement : 


Y+Dy=»+cL+c,M+c;Q 


Nous choisissons donc une formule à trois (quatre) approximations 
(voir J. Legras, p. 177) et nous déterminons les constantes a; b, c, de 
telle manière que : 


e l'écart entre la valeur exacte y et la valeur approchée soit infiniment 
petit d'ordre 3(4); 


e les trois premières approximations des deux formules d’ordre 3 et 4 
soient identiques. 


La quatrième O servira de première approximation L du prochain pas : 
Xo+h — x,+2h 
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Fehlberg propose le schéma suivant : 











a; b, C c; 
214 533 
L 0 0 — 
891 2 106 
1 1 1 
M : — — 0 
4 4 33 
o 27 — 189 729 650 800 
40 800 800 891 1053 
214 1 650 — 1 
O 1 — — — S — 
891 33 891 78 











23L M 350Q O 


BETETE TANT 








Celui-ci ne diffère guère d’un schéma proposé par Mineur et donne fort 
avantageusement l’erreur systématique E par pas. Cela nous permet de 
contrôler et d’ajuster la valeur de h au cours du calcul sans perte de 
temps. Si E est trop grand, |[E|>P, P erreur locale, on rejette la valeur 


h 
approchée y, + D y et on remplace h par 2 puis on recommence le calcul 


P 
à partir de y,. Si E est de très faible valeur, ES on remplace h par 


2h. 


PROGRAMMES 


19 REM ##*% Resolution d'un systeme d’ e 
quations differentielles d'ordre un *#* 
24 INPUT " Nombre d’equations n = ":H 

34 INPUT " Intervalle Ca. = "AE 

4@ IF A>=B THEH PRINT "aï=h7%": GOTO 34 
54 DIM YCMSUCHD, SCNS ICNDECNS,LCMNS, M 
CNY: OCND, GCHD 
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69 FOR K=1 TO H 
A READ IfK)D: REM ##* Valeurs initiales 


4#X 

84 YCKS=ICKO: SKI 

94 HEXT 

198 INPUT " Pas d’integration h = ";:H@ 


119 IF HO<=4 THEN PING: GOTO 144 

124 INPUT " Erreur lacale Pz "5 P 

139 IF P <= A THEN PING: GOTO 128 

140 *0 = A:H=H9:#=X4:G=04 

159 GOSUB A8 

164 FOR K=1TO N 

17A LCKI=H#UCK 3:ECK = 

1859 HNEXT 

199 PRINT" Valeurs %, 401),...:4Cn) :" 
208 PRINT #9: 

219 FOR K=1T0O H 

228 PRINT " ";:I<K5: 

238 NEXT 

24 PRINT 

250 IF x@>=B THEN END 

e6A FOR K=1 TO M 

278 SÉK) = SK) +214#LCK32/891:ECK3I=ECKD 
-234LCK 23/1782: YCK3=ICKOI+LCK 3/4 

288 NEXT 

234 #=x+H/4 

G9A GOSUR 744 

314 FOR Kz1 TO HW 

324 MKOISHKUC KO: SCKIESCKIHMCK I-AB:ECKDd= 
ECKO+MCK 1/33 

325 YCKIelICKI-ISSELCK 1-2 Q00+72S4MCK 5/20 
329 HNEXT 

344 #=#+M#17/49 

334 GOSUB 798 

360 FOR K=1 TO H 

379 GCKI=HKUCKOII SKI 
ECKOSECKD-3S54#0CK 0/11 
380 HEXT 

294 K=44+H#13/49 

408 GOSUB 744 

419 FOR K = 1 TO H 
429 OK = HXUCK): ECKÜ=E KO K 3-78 


K +6 SEEEC K 3-89 1 : 
YCKI=SCKD 


ss K 
Lea 
263: 


— 188 - 


£aa 
CL LOT 
814 


es 
EE" 


Fes bas pa pe pa pe 
+ ED M Bit 


ue 


fo 
Un 


30 
45 


HEXT 

RENM ##*# Frecision obtenue 7 ##% 
K=9 

KækK+1: IF K>N THEH 59G 

IF ABSCEGKODI2P THEN 694 

IF ABS ESK)3<P/144 THEH 469 
G=1:GOTO 469 

IF GzAG AND HCHO THEN HE=ZXH:GOTO 519 
IF H>=zHA THEN Gz1 

FOR Ke1 TO H 

ICK) = SéKD:LCK)D=C2-CoKOCKI:ECK ISA 
NEXT 

Geg: XAzx 

GOTO 2984 

HsH/2: #z#0: G=A 

FOR Kæ=1 TO H 

ECKY=@: SKI KO: LC KILL KI -E 

MNEXT 

GOTO 268 

REM #X# Systeme d’'equations *#%* 
Dei = CYCLE ZDD/C2#40 

UC23 = CVELI-VC2U CHA 

RETURN 

REM ### Valeurs initiales 4m(19.... 
2 #X*# 

DATA 1.9 


FEM ŒÆ#YX Trace d'une solution ##% 
PEN EE de l'equation differen-#&## 
FENM CHE tielle y''aefuuy ut TEKX 

REM ### uCaisun, u’/Caizum, acxéhE## 
FRIHTMytS œ =g jt 0" 

LPRIHNT Muyft = sin" 
FRIHT'Intersalle Cab = ":: INPUT A, 
PRINT "CA; ", "pm: "2" 

LPRINT "'luterwalle Ca,b9 = "i"c':Ai" 


;8;"2" 


IF A=E THEN PINS: GOTO 284 
PRINT "Yaleur initiale 4039 = ";i:IMF 


UT Y@: PRINT “4 
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45 LFRINT "Yaleur initiale sai = "; “4 


4 PRINT "Valeur initiale g'<aiïs ";:1INP 
UT UG: PRINT LA 
21 LPRINT "Valeur initiale s'faiz "; 4 


5 LPRIHT:PRIHT 
64 INPUT "Trace uk ou 4'CKIT 1/2 "AG 
: TeVALE HS) 
FA IF A$<>"1" THEH IF AS<S"2" THEN FING 
GOTO 64 
fa IF C$="0" THEH 85 
#4 CINy<Zz@ai 
85 M1=Yg: M2= Va: VC Ge YG:IF VALCASI=Z TH 
EN Mi=Ug:M2=U : 4/66 3=UG 
99 REM ###% Pas d'integration h=Chb-a3/28 
AXEX 
144 FEM ### xCkÿ= a+k#Chb-a3-244, VCk = 
éxCkIM Qu Yg'CMCHKOD KEK 
11 REM ### VCK2D Fichier ##*# 
124 H=CE-H31/-209: *A=H:K=4 
129 FRENM ### Methode de RUNGE-KUTTA ##%# 
149 “eg: Y=YAG:U SU 
154 GOSUE SAg 
16 Li=H#F 
170 X=#0+H/2:Y=YA+HEUA/Z II <UA4L1/2 
15 GOSUE 5A8 
199 L2=<H#F 
CAD Y=YA+HEUG E4+L1LHH/d: U sIA+L2/2 
218 GOSUB 54% 
224 L3=H#F 
239 #=#094+H: Y=YG+HEUR+MELZ-2:U <UG4+L3 
24n GOSUR SA 
254 L4=H*#F 
c6A <= 4+H:YOG="G+HELA+HECLI4+L24+L33/6 
274 UA=UR+#ELI+PÆLZ+TEL34+L 40/6 
284 K=K+1 
294 IF A$="2" THEN 339 
249 ne IF (A<M1 THEH M1=Y 
319 IF YA>MZ THEN M2=YA 
224 GOT 255 
234 YéKielG: IF UAR<Mi THEN Mi=Ua 
349 IF UG>M2 THEH M£zllf 
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tm ig D JM CAE 
DNA ERS RTE ET EN EURE ECr) 


à B BB 55 C0 QU (9 C9) CGJ CO 


45 
44f 
445 
45û 


IF K<290 THEN 14% 

IF Cé="Q" THEN 379 

HIRES: RENM ### Trace ##*# 
Ma=n1:Md4=Mz 

FOR Ke TO Z4a 

Mas LéRKk M4 VER 6 14-13 

CURSET K,W%:1 

HEXT 

FENM ##% Axes ###:1F C$="0" THEN 475 
Wa = 16A#CMd4-VC 045376 M4-1M2 5 

CURSET 4.V%,1: DEAU 705,4, 1 
CURSET 4.4,1 : CFAU 4.160,1 
CURSET 19, 179.@ 

HS = "Courbe gx) — trace M:IF T= 


£ THEH H$ = "Courhkhe u'fxi — trace " 


455 
44 
465 
4Fa 
474 
les 
473 


:":M1; 


476 
M2;" 
480 
Urs 
438 
EG 
495 


il 


DE Di M M MS MU LA 
LR Ride Tr iQmiS 


LS 


RO 2 UNS D io < ans 


FOR Kzi T0 LEMCNHS) 

CHAR ASCCMIDSCNE.K,1533,8,1 

CURMOY 14.4,8 

NEAT 

RE POKEZTE.SE Fermet d'ecrire dans 


3 lignes en bas de l'ecran Haute Res 


POKE 26: S6:PRINT "Yaleurs extremes 
LL " ; [ue ; LL LL i 

LPRINT "Valeurs esxtremes :"i:mi;" "ts 

M:LPRIHT 


PRINT"'Affichage du tableau des vale 


Q/H?°;::GET À$ 


IF A$= "H" THEN GOTO 644 ELSE GOTO 


RE ### Equation différentielle 
F = FAT UT ### 
=-S IMC 
RETLRH 
LPRIHT " Tableau des valeurs «Ck où 


Cyéhui" 


FOR K=Q TO 294 

LPRINT KA+KÉXCE-H1/2004, VCK5 

HET 

GOTO 664 

PRIHT" Autres valeurs initialest C0 


‘GET C$ 


IF C&="0" THEH 48 
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É6A PRIHT" Recopie ecran sur imprimante 
+ CO/Hi": GET DS 
FA IF Cf<>"H" THEH GOTO 1444 
2 EH 
194 RENM #X#X Imprimante XX% 
AA FRINT'AJOUTER SOUSFROGRAMME RECOPI 
E ECFAH HAUTE RESOLUTION" 

1495 MAIT 544 

1414 END 


De 


in 


UTILISATION DES PROGRAMMES 


Tracé d'une solution de l’équation y”=f(x, y, y’) 
On introduit la fonction en lignes 500 et 19 : 
500 F=f(X, Y, U) 
F=yX=x Y=yU=Y 


Le programme demande (lignes 20, 40, 50, 60) : 


e l'intervalle [a, b] (INPUT A, B); 
e les valeurs initiales y(a) et y’(a) (INPUT YO, INPUT UD); 
e de tracer y ou y'(1/2). 


(b— a) 

200 

Un fichier de 201 valeurs approchées de la fonction y ou y’ est dressé 
et les valeurs extrêmes M1 et M2 sont retenues. On peut choisir entre : 


Le pas d'intégration est fixe : h = 





e l'affichage du tableau des valeurs approchées (ligne 480); 


e le tracé sur écran (les axes passent par le point initial). 
Il est possible de superposer plusieurs courbes si les valeurs extrêmes 
des courbes suivantes ne dépassent pas celles de la première à tracer. 


Une mise en page sur l’imprimante est proposée. 
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Résolution d’un système 


Il faut introduire en lignes : 


20 le nombre d'équations N, 1=N 

30 l'intervalle [a, b] : A, B, A<B 

100 le pas d'intégration h : HO 0<h<0,2 
120 l’erreur locale p : P 0<p<107* 


les N fonctions du système 
700 UL=É YO. LL TON) EL UN = RE YEN ON) 
les N valeurs initiales y0 (1), …, yO(N) 


810 DATA... 


EXEMPLES NUMÉRIQUES 


Pendule pesant 





Coaurbe DR KE Où — trance 


y = sims 

Intervalle Ca.) = C6 14 1 
Valeur initiale 4623 = 
“aleur initiale s°'C4a2= 1 
Waleurs extremes :-1,094717477 


Walgur initiale 4ia4a1z 4 
Maleur initiale 4'éais ,S 


Waleurs extremes :-,545329743 
“Maleur initiale 4éai = 
Maleur initiale 4'6az —,75 
Maleurs gxtremes :-, Fé8re "ali 


Pont de Wien 


J 
LS 
AA] 
J 
ii 
EE 
& 
Li 
Eu 





J''+54"01-1,133/coshf2c y 13+yzg 
Intervalle Cab = C4 ,34 71 

Valeur initiale gai = ,4 

Valeur initiale s'faiz .2 

Valeurs extremes :-,T46256545 .746483 
CT RETTTEA 1 bg RUN Eu 

y'é Dee l ut Di CZ 

Hombre d'equations n = 7% 2 
Intervalle Cab] = 7 1,2 

Pas d’inteSration = FT ,45 

Erreur locale P= T7 1E-47 

WMaleurs %, 46117, 3402 

1 1 fi 

1.45 1.924995a4 .A24444174f 

Lai 1.H4996234 .AA4FEGIEZTE 

Lt 1.47485197 .HÉS9G9GSS851 

1.2 1.49975517 .n214142399 

1.25 1.1245A92E€ .112926155 

1.3 1.14919395 . 1219268223 

TL: 1.17378116 .151121975 

1.4 1.19526454 . 16727246 

1.45 1.22264431 . 15792763 

lie 1.24694413 . LASSLZLSAS 

LS 1.27145384 .2eréaerel 

1.6 1.299487 rÉ .239352953 

1.65 1.:1900486 . 2336024286 

1er 1.34220461 .2r1554529 

ir 1.2664969Z - 2HT1ES ZE 

1.5 1.39405245 22429397 

1.25 1.41354452 3172386982 
1.3a8aaaa1 1.42622208 . 2SCAEAALS 
1.95anaan1 1. 46005068 HAE 
2. 0444441 1.45315439 .:6h6185419 


Hombre d'equations "n = 
Intervalle Cab] = 7 1. 
Pas d'integration kh = 7 
Erreur locale P= 7 1E- GT 

Valeurs %s Clisson EN) à 


ê. 
1. 


F9 F2 1) 


1 1 q 

1.95 1.92439344 . 9244494174 
1.1 1.84936234 .AA476916235 

1.15 1.47435197 . 1699933831 
1.2 1.999723317 .4914142339 
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10 '= CLXYA 
OGT= x = 


ü= LY 
ART) AX 7) /.AD'TTI+A— AT 


NAN Tee 
LULU 


A#T—) SX7+T+( 
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16. ÉQUATIONS DIFFÉRENTIELLES 


AUX DÉRIVÉES PARTIELLES 








ÉQUATION PARABOLIQUE 


C'est une équation du type conduction de la chaleur : 


40. #0 


— = = 
OT AX" 


qui donne la température 6 en un point d’abscisse X, au temps T, d’une 
fine barre métallique qui ne transmet la chaleur que par ses extrémités. 


Il est intéressant dans un premier temps de normaliser cette équation 
de façon que toutes celles qui ont la même forme se résolvent par la 
même méthode. Posons : 


LT EE 


X,, et 6, étant des constantes pouvant avoir une signification physique. 
Alors : 


00 90 du 90 1 
= —X—=—Xx 
oX du OX du X, 





X, ou 





90. à (£ __à f1 20\ du 1 920 
äX? x el ) oX X2 ou? 
D'où, en reprenant l'équation de départ : 


à 1 9° 
2, F=ax xx 0, 


OT Xÿ du 
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aT Le 
Enfin, posant a on aboutit à : 
0 


pd 

ot ôu° 
Nous allons maintenant évaluer numériquement ces dérivées partielles. 
Par définition d’une dérivée : 

dp_p(+k)—g() 

ot k 


en appelant k un accroissement de t assez petit (nous précisons plus loin 
cette notion). 


De même, nous pouvons écrire : 


8p _p(u+h)—œp(u) 
ou h 


et tout aussi bien : 


dp_ p(u)—g(u—h) 
ou h 
Or, par définition : 


ap d@p 
vs (ee (u—h) 





au? h 
D'où : 


dp_p(u+th)—2p(u)+g(u—h) 
ou? h° 


L'’équation de départ nous conduit ainsi à la relation : 


pu t+k)—-p(ut) q(u+th1)—2p(u t1)+pu—ht) 
k u h° 
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_ 








(3 e e e e e 
, e e e e e 
k 
| u 
Figure 16.1. 
h 
œ(u, t+k) 
k 
œ(u—h,t) œ(u,t) p(u+h,t) 
k k [ < | 
[3] [2%] [Lr2 
| 
| | 
| | | 
Don oh | 
Figure 16.2. 


Considérons maintenant une maille repérée en abscisse par u et en 
ordonnée par t. En chaque point de cette maille, de coordonnées u = ih 
(i: entier positif ou nul), t = jk (k entier positif ou nul), existe une 
valeur de @—œ(u, t). 
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Or, l’équation précédente nous permet d’évaluer @(t+k) au point 
d’abscisse u, connaissant œ@(f) en u et @(u+h}), œ(u), g(u—h) à 
l'instant t. 

On peut ainsi écrire, de façon plus explicite : 


k 
PU t+k)=pu tr) + (putht)—2p(u t)+plu— ht) 


(Figure 16.2). 


L'évolution du système peut être ainsi maîtrisée instant après instant, 
connaissant : 


e les conditions initiales : répartition des températures le long de la 
barre à l'instant initial t—=0; 
e les effets de «bord» : températures imposées à chaque instant aux 


deux extrémités de la barre. 


Tout le problème est maintenant de réaliser un choix judicieux de h et 
de k. Les variables du système étant réduites, il paraît raisonnable de 
choisir h et k petits devant l’unité. Il a alors été démontré que la méthode 
n’est convergente que si : 


ce qui est fort gênant, parce qu’un h de l’ordre de 10? conduit à un k 
200 fois plus petit. 


Application 
La loi étudiée, de la forme : 


äp 3° 
ot du? 





est très généralement connue par les physiciens sous le nom de loi de 
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Diffusion (ou seconde loi de Fick); elle exprime la conservation du 
nombre de particules; elle s'écrit : 

on don 

dt du’ 
dans le cas unidimensionnel : 











Figure 16.3. 


Remarque 


Une solution avec n(0o, t)=0 et n(u, 0)=1 est n(u,t) : = 


D 
* (V4) 
21 
= fonction d'erreur (voir le Chapitre 5). 


Soit n la concentration des molécules d’un gaz dans un conduit 
parallélépipédique. Si cette concentration n’est pas homogène à 1—0, 
elle est appelée à varier en fonction du temps; cette évolution est régie 
par la loi précédente. 

Considérons l'exemple suivant : 


Soit un conduit rectiligne reliant deux récipients entre eux : 








ni no de 





| 
| 
+++ —+ u 


AE 





Figure 16.4. 


A t=0, la densité de molécules est homogène dans le tube et vaut n,. 
Elle est imposée égale à nr, (à gauche) et à n, (à droite), quel que soit le 
temps t. 
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Le problème est de savoir comment va varier, à partir de t—0, la 
densité n(u, t) en chaque point du conduit. Le régime permanent peut 
être connu par un raisonnement simple : 

Si, au bout d’un temps infini, s'établit un régime stationnaire, c’est que 
l’on a alors : 


on _ 
ot 


On en déduit : 


dn 


du? 
La densité va varier de façon linéaire à l’intérieur du tube de n, à n.. 


Méthode de Crank-Nicholson 


La méthode étudiée précédemment est qualifiée d’explicite dans la 
mesure où l'équation obtenue, permettant la résolution numérique, 
fournissait la valeur de @ en un point en fonction des valeurs déjà 
connues de @. Le problème pouvait ainsi être résolu de proche en proche. 

Cependant, la condition de convergence : 


impliquait d’utiliser un pas temporel très petit si l’on voulait obtenir une 
précision spatiale correcte. 

Crank et Nicholson proposèrent une méthode dite implicite, car elle 
faisait intervenir une équation de base ne fournissant pas aussi directe- 
ment la valeur de @ en chaque point. L’équation de base étant toujours : 


dp 9° 
dt du? 





> 


la méthode diffère par le choix de l’expression utilisée pour à 
u 
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On peut en effet décider que cette dérivée est la moyenne arithmétique 
des valeurs numériques qu’elle a aux temps ft et t+k: 


2 @p._ 1 
du? 2h? 


+p(u+tht)—-2çq{(u t)+œ(u—ht)] 








[Q(u+ht+k)—2q(u t+k)+q(u—ht+k) 


k 
D'où, pour l’équation générale (en posant toujours T3 et en ayant 


, dp plu t+k)—œ(u, 2) 
LOUJOUTS — = ———— |): 
ot k 


—Tp(u—ht+k)+(2+27T)qp(u t+k)—Tp(u+h t1+k)= 
T@(u—ht)+(2—27T)@(u t)+Tq(u+h,t) 


Il est alors commode de reconsidérer le réseau maillé défini précédem- 
ment. 


Prenons t—0. L’équation précédente, appliquée aux points d’abscisse 
u—=h u=2h,.., u—=nh—1, en appelant nh la longueur horizontale 
du réseau, fournit les valeurs de @ au temps 1t=0+Kk, en fonction de 
celles qui existent à 1 —0 (qui sont connues, puisque ce sont les valeurs 
initiales). Nous pouvons ainsi écrire n — 1 équations, correspondant aux 
(n — 1) valeurs cherchées de @ sur l’axe Ou au temps t = k. Ce problème 
peut alors être résolu classiquement (système de (n —1) équations à 
(n — 1) inconnues, voir le Chapitre 12). 

Dans une seconde étape, connaissant toutes les valeurs de @ à l’instant 
t—= k, il est possible de recommencer l'opération pour en déduire les 
nouvelles valeurs de @ à l'instant t=k<+k, =2 k, et ainsi de suite. 

Cette méthode, qui semble apparemment plus longue que la précé- 


dente, a l’avantage de converger quelle que soit la valeur de = 
Cependant, la précision de la méthode est liée à ce choix de 7. Une 
valeur telle que T = quelques unités peut convenir. 
En général, on choisit T—1, car cela permet d’annuler certains 
coefficients de la relation établie plus haut. 
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Là INPUT HA. H1. HE 
24 IHPUT'HONBRE DE POINTS" iH 

3 INPUT'ÉRELISION CL OHRER, SoeE 
bd Mel-H:REM PAS SPATIAL 

25 KeftMoe:REM FAS TEMPOREL 

6 LM H 5: DIE H 3 











27 HIRES 
HA 


dà FÜR J=1 TO 3 

SG HE 1 età 

Hä CURSET 1,193-Ac13.1 

EM MEXT I 

F4 HE eh] : A LH eHE 

CURSET ,199-f6 60, LCURSET LAN, LAS LE 1 





1 


Aù EUR H = & TO 1MGYSTEF 24 

5 FÜR UeH TO H+135 

GE FÜR I=1 TÙ 2% 

LE BE I ORAC LE KE HUB MRC HO LL BAC LAC [1 05 
34 MEST 1 


154 FÜR si TO 83 
168 HET 2=E T0 
1F4 HEST I 


163 HET 
154 FÜR I=1iT09a 
9. LURSET I,199-61 13:1 
156 HEXTI 
156 FRIHTH 
134 HEATH 


Programme BASIC 


LA LHELIT HE HI OH 

CUT LEA LA 5e D TES 1 AA 
ét LHELIT HE 

cé MIRE 

td 

HA RCI Lei TI 2 
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[tdi 

An CURE LE Li Lite Ta, 
4 HE 
CET HE 
a CUIRSET 

ral 

4 FOR Hz A TO LAGASTEF HA 
a FOR LE TO H4+134 

Lit FR Let Tir 

LEA OMC T et Le EE MO I CEE ME T4 et Æ AE L'HAE Led 
LS HÉAT TI 

1 ti 

LS FOR [el TE 34 

LE ET ieEe Li 

LA HEMT T 






LA, L'A-Aé Lti 2, À 


LE HEAT I 
154 FR I1=17033 
Fe ILURSE TE, 199-B60 2; 





156 HEATI 
159 FRIHTH 
14 HEXSTH 


PASSES 
a 


=! 
“ # 


& 





D 
LATE 
‘ 
LAC 
ue 
re 
* 
n 


pd 
A 
% 
4 
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Programme Pascal 





PR 
EE 





3 

DE 
DE 
ER 
DE 


S 
= 
— | 
L ; 
ua | 


in 
Le] 
Et 
nel 
pu] 
1 
ti 





z 
Lu] 
FI 
Fa] 
jf 


tt 

Es 

EN fe CU TL 
# 

fa ba 


{+ 
PROCEDURE 


SESIt} 
HRITE'C 
& IH se. 
:SABE, 


3 
T1 


mm 
[ra 
qu 





[TI 








Ne Pi a CN 





#a1 FE, fc 
SE, #F3i: 
"ENS; 
{+ 3: 
BEGIN 
BEGIN 


Ces CTI CT ee Ci PA ee ee ET be be 
x] ue =. 








$C es 





READ ik) 


RO LOC LAC LA LOC LAC LCI PEN ET ES ES EC 





CD F0 C0 CO CN 








“| 
DISTES 
SES LS) 
L.4 | Li 
[+4 (0 
Ita fl 
Hume 
À va + ra = 
LU A ee A 
IL Là LA 4H LE LA 


6e na EE DUT 





= 80 
N2=20 


N1=140 
NO 


2 


EQUATION HYPERBOLIQUE 


pas spatial 0,01 
pas temporel 10—* 


équation du type : 


C'est une 





tion de propagation (nous la supposons ici 


équa 


Il s’agit d’une 


normalisée). 
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En utilisant la même relation que précédemment pour exprimer 
numériquement ces dérivées partielles (pas spatial = h, pas temporel = 
k), nous pouvons écrire : 


qlu+ht)—-2@(u,t)+q(u—ht) q(ut+k)—-2p(u t)+q(u, t—k) 
h? =: k°? 


ce qui devient : 


pu t+k)=@ q(u—ht)+2(1—-æœ)q{(u t)+a piu+ht)—qpiu t—k) 


en posant : 
k 
a=— 
h 


Cette équation peut être représentée dans un plan paramétré en uett 
(voir la Figure 16.5). 

Il a été démontré que l’emploi de cette méthode ne conduit à une 
solution convergente que si a<1. 


t 








Figure 16.5. 
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Il existe une seconde méthode permettant de résoudre ce genre 
d'équations, en ramenant leur résolution à celle (simultanée) de deux 
équations différentielles du premier ordre. 


Il suffit pour cela de poser : 


ù ù 
= f(u t) et = g(u t) 


La relation : 








dp_dp 

dudt dtdu 
jointe à : 

ap à 


ou” ot” 








conduit alors à : 


df(u.t)_ dg(u,t) 


ou ot 
df(u,t)_ dg(u t) 
ôt  ôu 


Ce système d'équations différentielles peut alors se résoudre par un 
procédé classique. 


Revenons à la première méthode et tentons de l’exploiter. Il est évident 
que pour commencer un calcul, il est nécessaire de connaître la valeur de 
la fonction @ en quatre points (valeurs initiales) pour en déduire la 
cinquième (et nouvelle) valeur de @. 

Considérons alors la Figure 16.6. Elle pourrait correspondre à une 
corde fixe à une extrémité (à droite) et excitée à gauche par un vibreur 
imposant une oscillation verticale sinusoïdale. Il s’agit d’un problème 
classique dans nos classes de terminale, connu sous le nom de cordes 


— 210 - 


vibrantes. Un résultat classique est le suivant : si la longueur de la corde 
est égale à un nombre entier de demi-longueurs d'onde, alors celle-ci est 
le siège d’une résonance, résultat de la superposition des ondes directes et 
réfléchies. Il est alors paradoxal de constater que, dans cette situation, la 
source est un nœud et l’amplitude de la vibration aux ventres est très 
grande. La seule différence entre ce problème réel et la simulation que 
nous proposons est que cette amplitude aux ventres est certainement 
limitée en pratique par l’extensibilité de la corde. Cependant nous 
constaterons que, dans notre simulation, cette amplitude ne devient pas 
non plus infinie, cela résultant du conflit entre la source et le nœud 
qu'elle devrait être! 


Revenons à la Figure 16.6. Nous supposons la position des différents 
points (il y en a 50) contenus dans trois matrices : B(50) à l'instant 
t= —2k, A(50) à l'instant 1 = — k, C(50) à l'instant t=0; la corde est 
ainsi supposée immobile pour 10. A partir de 1 —0, ces données nous 
permettent de déduire les 50 positions des différents points, sachant que 
C (0) est l’excitation sinusoïdale (la source), et C (50) = 0 quel que soit t, 
puisqu'il s’agit de l’extrémité fixe de la corde. 























+ Points connus + Points calculés 
à l'instant 1 —=0 à l'instant 1 = k 
t 
| 
| 
introduction 
de C(0) 
1=k 4— ———— 
0 0 
à=0 dpt ++ —c(1) 4e —e—e + — AI) 
0,0 ,0 0,0 J0 
En Ok PE == +++ A (1) +—e e + 9 © e—4#-—B(I) 
0 10 
1=0— 2 & 40-40 40 40 0 ,0 40. B(1) 0 1 2 3 48 49 50 
1=0 1 2 3 48 49 50 
Figure 16.6. 
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Cela étant fait, nous pouvons passer à l'instant {= k. 
Les deux positions «anciennes » dont nous avons besoin pour effectuer 
les calculs sont par définition A (50) et B (50). 


Or, ce sont justement les matrices C (50) et A (50) évaluées à l’instant 
t=0 et (= — k. Il suffit d'effectuer le transfert des coefficients des unes 
vers les autres pour pouvoir effectuer une nouvelle itération, d’où le 
programme : 


1 INFUT"'PRECISION CHA 15 = "A 
2 LET H=1/54:K=A#H 

3 HIRES 

19 DIMAS 59: DIMES SA 1: DIMCE SA 
24 FOR 1=Q TO S4:ECTI=G:MEXT 

34 FÜR I=1 TO SA: ALI SED:HEXT 

44 Ce @i=S INC 2#PIET 3: CC SG = 

59 FOR I=1 TO 43 

69 CET O=ALAEAS IL IEEE 1 -AÆA TÆEAC TI + HEHEA 
CIt+12-B6 10 

FA CURSETI, 144-SS#Cc I 3.1 

S4 MEXTI 

98 FOR I=4 TO = 

199 BCI2=AC I) 

114 HCII=CC là 

129 MEXTI 

154 T=T+K 

144 HIRES:G0TO4f 


Il ne faut pas oublier, comme pour tous les problèmes de ce type, qu'il 
y à au départ un régime transitoire qui, suivant les conditions initiales, 
peut être très long. 
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PEU RALHEHUOES 
, 





@) 0,36<1<0,72 


ÉDCLEUL EUDDE CRC ESUT CE 






nee D 
ALLAN LA LE CPE PERS 


@) 1, 0,72<1<1,08 





DORA ENT HE ENEENE 
@ 22 
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0,5 
57 
0,01 


Il 


F2? 
| 








96 itérations) 
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è 
st É , 


. . . 
‘ 
COPPPS CEE CITE CIT ELLE 
, 








LOLLELT EE LESELE, 
LEPPCIT TPE Û 


4. 














D 
CR 


ne ee 


= 
[ex 
(a 
@) 
b 
Q0 
O0 
= 
où 
= 
[NS] 
2. 
© 
D 
a 
sr 





‘ 
Cr 
. 
nm 
. ; 
COPPPES LEP Cr . 
‘ 
... 
Un Gites 
LL 
AIR 
D 
AUILLLTEE 
. D 
DE 
M. 
‘ 
LE 
‘4 
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î 
ions COELETELEL ILE ELEC CLEPCETEETECEEE EEE EL EEE 


oO) 


H, 
Missoouss 


COTÉES TEE EEE EL EEE EEE CEE CEE E LEE 


fn. COLLI TETE EEE EE TEE disseiruréene ssesuse 


©) 


lu, 
d 


fn. [ETES CLLT EPP EEE EEE EE) doontoososne QE 


@ 
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COLUTTR 
“it. 


è 
“lite 
M 
‘tt 
: 
Hs CODE EPPEC ET EEE TELE TETE EEE 


HAUTE 
rl M, 
il 


LM hrs Booocosoooooosooss 
MUQUTE 
RL 
RL —— 


Fe 
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26 PROGRAM HY 
2 LABEL 3e ; 
CONST 








DNÉHEMAIETE. 


& I 
} à LAR I: TEGE 
RCE7 & KT: INTE 
RDE7 SE T,F:RERAL: 
ACEF iaë R:'ARRRBT IS, HIDCF RE 
AL ; 
RACS7 116 B:ARRBE IS, HIER BE 
PE 
LE 3» 
BDE7F i£Sg LC: ARABY [E,. RE 





FROCEDURE CCF"; 

BEGIN 

INLINE HFC.HREL,: 
HOI,HCE, HCD RAC ,+ 
L,RSE.HFai; 
END ; 
FROCEBURE 


Lans 


LD 
HUIT 
TE ON EN 

CE TD 
HH 
CIS te 





CL LE LEE 


ba pa la 


LT RTR 





BEGIN 
URITE: 


EX 
# 


3 6 Ga 


DODURIE 





Ti 

mm 

re 

I pe 19 K70 19 
Fa hs hi 








“# 
Le 
H 


ES rTrT 


H4 
Il 

Le] 
4 
CN] 











I: =8 
109: C [0] 
C [N] : =@; 
FOR I:=1 
AF BEGIN£SC 
AF CIIT:=iE 
il IJi+iiEs 
EREL IF T 
H à ; 
E 


PLUIE 


F 
a 
Et? 
B1i4s 2 
THEW FLOT 
5134 ë 
B1355 à 
=B LI] 
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1 


C0 EU CD CTI EN 11 € 





ïi 





em 
es 
Le ii 
MS: 
= 
27FS 
sers 
pen) rs 
==; 
De. Los 
RARES 
TL 
Er 
a 
= 








ÉQUATION ELLIPTIQUE 


C'est une équation de la forme : 


9° HE 
PE 


SE are ant À 
ou” y 


connue en physique sous le nom d’équation de Poisson. Plus générale- 
ment, on peut l'écrire : 


9 9° 
LR Tete y) 
ou 9y 


En utilisant la même relation que précédemment pour exprimer 
numériquement ces dérivées partielles (nous supposons la maille carrée, 


c'est-à-dire que le pas sur l'axe Ou, soit h, est égal au pas sur l’axe Oy 
(voir la Figure 16.7)) : 


O=qu+h. y)—2çq{u, y)+q(u—h—y)+çq(u y+h)—2ç(u, y) 
+œ(u, y—h) 


D'où : 


4qlu. y)—qiu+h y)—qi(u—h, y)—çq(u y+h)—œlu y—h)=0 


Considérons par exemple une plaque rectangulaire, pour laquelle la 
température des bords est imposée par un certain dispositif et sur laquelle 
nous disposons un certain nombre de points de mesure équidistants. 
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100° 100°  100°  100° 











11000: 100° 
| 100°® ) 100° 
100° D ) 100° 
0° 0° 0° 0° 
*<—— 3 + ———>D<{ >< << +» 
h h h h h 
Figure 16.7. 











+ uU 





Figure 16.8. 
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La relation précédente impose une relation entre certains groupements 
de points. On peut symboliser cela par le dessin suivant (Figure 16.8), qui 
est une représentation spatiale de cette relation. A chaque point 
d’abscisse u et d’ordonnée y est associée une valeur de la fonction 
p(u, y), elle-même reliée aux valeurs de @ pour les points entourant 
celui qui est considéré. 

Revenons à la plaque rectangulaire définie précédemment. Supposons 
que tous les bords soient maintenus à la température de 100 °C, sauf le 
bord inférieur qui, lui, est porté à 0 °C. Intéressons-nous à la température 
de chaque point intérieur (!). 





Figure 16.9. 


Considérons le point (1). La valeur de la fonction @ en ce point, soit 
®,, est telle que : 


4p—1X(100)—1X(p;)—1X(p)—1X(0)=0 





1. D. Herrmann : 40 BASIC Programme, Vieweg, Wiesbaden, 1983. 


G. D. Smith : Numerical Solution of Partial Differential Equations, Oxford University 
Press. 
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soit, en ordonnant : 


4ap—I1Xp +OXp,+OX p—I1X p+..—=1X(100)+1X0 
=1 X(100+0) 
Continuons; pour le point (2), nous avons : 
4 D — Pi — Ps — P3 —0=0 
D'où : 


—1Xp+A4AX p—1Xp,..—1X p,+..—=1X (0) 


On conçoit donc qu’à chaque point intérieur sont associées : 


e une série de coefficients; 


e une température de «bord» égale à la plus proche considérée, sauf 


pour les angles ((1), (4), (9), (12)), auxquels on associe la somme des 
températures des bords de coins les plus proches. 


Ainsi : 


u, — 0+100 
u, —+ 0 
u; —+ 0 
uy — 0+100 


us —> 100 


us ——+> 100 +100 


no d00 Æ 100 
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On peut ainsi écrire, en tenant compte du fait qu’à chaque point du 
plan est associée une ligne de coefficients : 


100 4 — | 0 0 — 1 0 0 0 0 0 0 0 gi 


0 0 — | 4 — 1 0 0 —1 0 0 0 0 0 g 
100 0 0 — | 4 =] 0 0 — 1 0 0 0 0 Pa 


200 0 0 0 0 —1 0 0 —1 4 —1 0 0 Pu 


0 0 0 0 0 0 —1 0 0 —1 4 —1 0 Pro 
0 0 0 0 0 0 0 — 1 0 0 =] 4 — 1 Pi 


200 0 0 0 0 0 0 0 —1 0 0 —1 4 Pi: 


Le problème revient donc à trouver la solution de ce système 
d'équations. 

Il existe alors deux types de méthodes applicables à la résolution d’un 
tel système : 


e la méthode directe (résolution du système d’équations, qui peut être 
très longue; voir le chapitre concerné); 


e les méthodes itératives qui, partant d’un ensemble de valeurs @, assez 
quelconques, conduisent après plusieurs passages à une convergence 
éventuelle vers la solution. 


Nous décrivons ici deux de ces méthodes : 


Liebmann non extrapolée ou méthode de Gauss-Seidel 
(méthode de relaxation) 


Considérons le (n + 1)*"* passage; nous pouvons écrire : 
re [ (n4+1) (n) + (+ E (n) 
Pi. “A Pi-1.; PitiT Pi ji Pi j+1 
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Cette équation permet d'évaluer au (n +1)“ passage @,, connais- 
sant : 


e Piyi jet Pi ;+1 évalués au passage n précédent (points cerclés sur le 
dessin); 


e Pij-1 et Pi-1, évalués pendant ce même passage n +1 (points 
entourés d’un carré). (Voir la Figure 16.9.) 


g., est aussi estimée à l’aide de valeurs de @ les plus récemment 
connues. 


Il a été montré que cette méthode converge. 
On peut évaluer la précision à la (n + 1)*" itération, en calculant le 
résidu R, ; au point (à j) : 


— /n(n+1) nn) (n+1) (n) = On) 
R;= pi; + pi, ; + Pi j—1 + péri 4 pi; 


quantité qui est évidemment nulle pour /a solution. 


























Figure 16.10. 
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Exemple 


Considérons une fonction @ à deux variables, qui obéit à l’équation : 


dp à 


du? y 


= 0 








Dans un réseau rectangulaire (voir la Figure 16.10), cette fonction 
prend les valeurs : 50 pour les bords gauche et droit, 100 pour le bord 
supérieur, 10 pour le bord inférieur. 

La matrice P (I, J) contient à chaque instant les valeurs de la fonction : 


pri (1, J) 


Q(I, J) contient les valeurs antérieures @("? (I, J). A la fin de chaque 
itération, Q est remplie avec les éléments de P. 

Les résultats (valeur de @ et des résidus) sont imprimés tous les cinq 
passages. 

Pour commencer le calcul, il faut donner des valeurs arbitraires initiales 
aux éléments de Q. Ici, nous avons pris Q(I, J)= 50, ce qui semble 
raisonnable. On aurait pu également faire un calcul préalable, en 
n’évaluant la fonction @ que pour un point sur deux (16 points à calculer 
au lieu de 81). 


à CIN FC19,163%:DIM 614,145): DIM RC19.1 
A1: H=G 

DA REMESSrneennes smes 
cA RENM CONDITIONS AUX LIMITES 
DD PEMS=emsunmnss nn sers 
44 FOR I=4 TO 14 

ag FCI.G3=14 

a QCI.Gr=14 

SA MEXT 

FA FOR I=4 TO 14 

64 FCI,142=148 

Q4 GET. 141-144 

9: MENT 

144 FOR J = @ TO 14 

114 FCG..11=5 a 
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Gé 4, 11254 


115 “12 

124 HEXT 

134 FOR Je: TO 14 

14 FC1A..13=5 A 

SA QC14,.11=5 4 

168 HEXT 

LA REMeee ne ne cs ce ee ne ace nan ee ee 
21 RENM ETAT IHNITIAL QUELCONQUE 
ERA REMeee nan ce ne nee ae ae ee ae ae ae a 
224 FOR I=1 T0 3 


Ed4 FOR Jel 
ss GI, Je 
264 HEXT 
274 HEXT 
BQE REM-—--———- 
319 REN CALCULE 


TO 3 
ft 


229 FOR J=1 TO 4 

Es FOR Iz1 T0 3 

SA PCT. Joe, PK PC I-1. +0 T+1. JP. .J 
AH T IA M 

HEAT TI 

MEXT .] 

H=t+1 

FRIHTH 


ire 18 20 0 TN I 
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PEN EVALUATION DU RESIDU 
ABQ REM— 

434 Ie INTE FHDC 1 3941 

448 JG=IHTÉ RHD 1 VE] 

4SR R=Pé TA-1. JA QE LL, GARCIA A1 34 
QC IG, JE+1 1-0 IQ, JE 

460 FEsR Pi 16, 0 

47 PRINTRS 


ai FEM ITERATION 
Art EM 
SA FOR Izi T0 
548 FOR ei TO 
4 RET, MiePe T1. Joie T1, PCT. 1 +0 
TU dE TT, 15 

SA DT. J'ieFé TT, lt 


ni 


— 233 - 


Fa 


SAGE NUMEr 


LUN 
a 


MEXT I 
HEAT TI 


Mini 


=, 
En 


_GOTE AE 
FOR .I= 
FO I=A TA 
LFETNT TAF 
HET I 
LFEIHT 
HET I 


DR ét db pb be pes pt is pe EN TEA EN CN 
He a He pe pa pee pes ps LÛR LA Of T3 


D EN BR 9 fi me D Sn à 
CSSS SAS AE 


RETURH 





Œ 
ct 


1 


THE: IHTÉRET. JI#14 3-18; 


LPFRETAT:LPETHT" Pa: 
LFRIHT:LPRINT:LP 





= 934 


(gl 
14 TO G STEP -1 


44 
T: 


PRIH 


T 


Rx 


E 


IF H=1 THEH GOSUE 11644 
IF INT HAS 


Muret 


LPEINT 


mi 


CA AE 


1 


Rumis 


_.i 


= Lo ee ei 


NS 


it 


Lu 


IT 


H-S THEH GOSUE 114% 


";H 


Tin 


ji il 


Ti 


AS SRE D ic 
Te TS D ET EE 


AARRT 


DE 


Sn Et CO GO DR OU D D D EE RE EE D SE 
FEU to 4 D I M tt ei mt ei mt m4 
x PRE SI es de en Dose sms 8 QU ess GRR DDR ss DRE 


DUT Fe Fe UT En NES 


us en 8 es 0 






eu RSS DRE 


Fi 









a EE UT UT D OU 4 Par 


Das se GS UD ss ss ss 





ni 
5 
1 
1 
ï 








Dam D + Qi AO M <f EN 4 


sn Hess... sw inst 





D EU 4 
GS 





M es + + Aa A Fe HE du 


Ex) D eme 8 = Gers To La a Us de a 









EU UT LE Fe QD OU 



















Bo bal Rss rsss se .& 
} 
L mi brx.| 
+ EX] St sum 8 si 
4 
ba) ns DA] 
_ ee Re} 
mi + 
UD © o Mu us ss [nl Du RE 
L L L 
re ü En 
= = mmMNRiU en mt = 4 et ei 
Ex osseuses DE es Um EU 
& l ÿ à 1 ü 
œ m 
; : 
ut 
TL = ei ee —. ne ui ee 
Re D En CR EU SE Et EE EU RE A EEE EEE ET) es SRE ARE 





— 235 - 








ñ Fi 4 E] A ñ (si) 
[E F [a F4 [ns] ni] A 
f A ñ (al (a f ñ 
ñ fi “1 ra [a] [a] £ 
Ca [a] il sl 4 (a) La) 
F1 (a) il 2 sil ] f 
[al [a] 5 ue “À [a] f 
A (E] sd s1 SL [a] A 
fi a [a 1 . al .l Fi ñ [s! 
(] F4 (El A A [e] ñ A f ñ [x] 
(e] A A [e 4 A 4 [s] fi 4 (e] 
Passañe numero M 

[) Fi Hi 1 ñ ä ä [a] [al [a] ñ 
fi F1 F1 Fi [a] [a] [a] [al ä (s E 
F1 [el [a li A E] (S. A (| 4 (] 
F1 [) (a) 1 .] ml ‘i [a] [) fä f 
(] A Î 2 À HA .1 3 [a] A La 
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A [a] a (a) (a) [a] f [a] [a] (a) A 
[a] [el (a) (el [al [al [al g [a] La] La) 
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ñ (a) 4 [al ñ A ä [al (a) [a] A] 
[a] [a] a [a] ä El [a] A (a) g [z 
4 [a] [e] (E) [a] [a] (El (a) a Il f 
ñ (e] [a] [a] a [ul] 4 G [a] (al ä 
ñ [a] La La] [al a [a] [a [al f g 
f A (s) A a] [a] a] g Li] f (s] 


FassaQe numero 4 


Remarque 


Ce changement de la ligne 1120 permet d’imprimer les valeurs de : 


11524 LPRIHT TARCTÆR CS TNTERET. NET 14: 
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min 









1 
CRE 
3 €: 
TS 3 F2 
sà € 5 1e 
SE F 8 4 SE 
la 9 .4 4 4 
pal D .Ë 4 À 
54 «5 7 & 3 
DA .H dé sl 3 
54 14 14 1 
Passage numero df 


Liebmann extrapolée où S.OR. (successive over relaxation) 


L'itération précédente peut aussi s’écrire : 
(n+1)— pin) 1 (n+1) (n) (n+1) (n) 4 ot") 
Pi; = Er + LOT + Pi, + EE + Pin — 4 pi 


Mais le terme entre crochets est précisément ce que nous avons appelé 
précédemment R,; ; résidu au point à, j : 


1 
pie pi + Ru; 
Dans cette méthode, on est amené à poser : 
(n+1)— pin) 1 
pi =; + 1° R;,; 


1 1 
où G «) est un terme positif compris entre à et 3 en pratique. 


Remarquons que la méthode précédente (Gauss-Seidel) revient à prendre 
=]. 


1 
En posant nr &, il a été montré que la valeur optimale de a est la 


plus petite racine de l’équation : 
az? —4a+1—=0 
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dans laquelle : 


T T 
Z = COS — + cos — 
Im JM 


iu = valeur extrême (maximale) de l’indice i; 


ju = valeur extrême (maximale) de l'indice j. 
Exemple 


Reprenons l'exercice précédent et résolvons-le avec cette nouvelle 
méthode. Nous avons : 


iu= ju = 10 


L'équation du second degré précédente a deux racines 0,73 et 0,38. 
1 
Cette dernière sera notre coefficient à &. 


Le programme, modifié en conséquence, conduit beaucoup plus 
rapidement à la solution. 


5 DIM P619.199:DIM QC14,1@5: DIM R<10,1 


B:N=9 

19 REM--mmma soon oo mt ne 
24 REM CONDITIONS AUX LIMITES 
DO REM-—-- une ne not ee mme moe 


49 FOR 1=9 TO 19 
59 PCI. :=1n 

55 AMI.a1=18 

6A MEXT 

79 FOR 1=4 TO 19 
8a PCI,182=164 

94 Q1.143=189 
93 NEXT 

194 FOR J = 9 TO 18 
119 P6Q,.J)=54 

115 960,.13=59 

129 MEXT 

139 FOR .J:9 TO 18 
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149 PC10,J5=52 
159 QC 14..13=59 
16A MEXT 


239 FOR I=1 TO 3 
24 FOR J=1 TO 3 
259 QCI,J2=59 
264 NEXT 


M 
314 REM CALCULS 

329 REM------ 

330 FOR J=1 T0 9 

349 FOR I=1 TO 3 

358 PCI, Je I, JHE#E PE I-1,.J 04 I41, Je 
PCI, J-1 +00 I, J+13-4#QC 1,133 

369 NEXT I 

379 MEXT 

389 N=N+1 

394 PRINTN 

490 REM-————— 

419 REM EYALUATION DU RESIDU 

420 REM-—————— 

436 1G=IHTE END 1 24941) 

448 JO=INTE RND 1 2941 1 

459 RaePe1@-1,.J2 +06 1041, J@ 3+PC IG, JA 1 + 
QC 14. J9+1 3-4#0C IA, JA 

464 RS=R/PC 10. JE 

474 PRIHTRS 

SAQ REM-—------- 

“19 REM ITERATIOH 

G2Q REM------- 

529 FOR I=1 TO 9 

54@ FOR J=1 TO 9 

SdS RCI: YePC T1, J40C +1, IPC I. 1 340 
+1 3-4# 06 1, 0 

55Q OC, J>=P£ I, 

562 HEXT J 

579 HEXT I 


nt 


Te 


SEA IF H=1 THEH GOSUE 114 

SES IF INT éH/SeH-S THEN GÛSUE 1144 
534 GOTO 3984 

1104 FOR JI=1G TO G STEP -1 

1114 FOR I=A TO 1 

1124 LPRINT TRECI#E M: IMTORET. Jik14-14:; 
1134 HEXT I 

1144 LFRIHT 

1154 HEZAT I 

1164 LPRIHT:LPRINT'"'Passage numero ":H 
1165 LPRINT:LFAIMT:LPERINT:LPRIHT 

1174 FETLUFH 


MERSION COULEUR 
9 DIM P611,119:DIM 611,119: DIM RC11,1 


1 5: H=û 
1Q RE 
24 REM COHDITIOHNS AUX LIMITES 
BA REM 
44 FOR I=9 TO 11 
SA PCl.43=10 
55 QC 1.G)=14 
6A HEXT 
74 FOR I=f TO 11 
EM Pl, 113=180 
94 ECI,117=180 
95 MEXT 
aa FOR J = A TO 11 
114 F6, .J3=5a 
115 069,125 
124 HEXT 
154 FOR Jeë TO 11 
149 P<11,.13=50 
154 QC11,1=50 
168 HEXT 
BGQ REM 
214 FEM ETAT IHITIAL QUELCONQUE 
DEA REM-— ne 
234 FOR I=1 TO 16 
244 FOR J=1 TO 14 


AR = 


54 QT,Ji=Ssm 
A NEAT 
Pa HEXT 


DGA REM-----—- 

318 REM CALCULE 

SEA REM------— 

334 FOR Jei TO 18 

34@ FOR 1=1 TO 18 

359 PCI. Je, 2SKC PE IA, 9406 141, J4PC I, 
1 +0 413) 

364 MEXT I 

374 HEXT J 

580 Heh+1 

AUD REM-n nn 
418 REM EVALUATIOH OÙ RESIDU 
ABG RM 


434 IG=IHTERHOE 1 141041 3 

440 JG=THTE HD 1 41041 5 

454 RePCIG-1, IQ 406 TD+1,.1@ +6 T4, 4-1 + 
QC TG, l+1 -d£cé TG,.14 5 

HS=R/PC TG, IG 

RÉF 

REM ITEFATIOH 


ninin R& 
\ù te I M 
oi 1 


234 FOR I=1 TO 14 

944 FÜR Je1 TO 14 

945 Ro, Jo=Pé I-1, 4 GC T +1, JP TI, J-1 140 
ÊT. J+1-dÆ0ic TT. 15 

aa4 QT. IsPT. 3 

364 HEXT .] 

374 HEAT I 

Di IF H=1 THEH GOSUE 1149 

5 IF INT CMS eH/5 THEN GOSUR 1144 

14 GOTO SAR 

FOR I=4 TO 11 

FOR KG TO à 

CURSET 124+I#1S+K#€6, 14,35 

FÜR J=11 T0 6 STEP -1 

Wi=C PCT. 11-1447, 99/94: M2=U1-INTCN 
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Ce 
nn 


Bones 
Pres be iQ iT 
Ron D En S 
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DAS 


1125 M=16+IMTCUI 
1159 IF CKeli AND CUZ%,.665 ANDCU<232THE 
H Hesl+1 


11525 IF CKet5 AND CH2<,2353 AHOCH@3 THE 
H W=h-1 

114% FILL 15.:1,4 

1145 HEXT 

1159 HEXT 

1155 HEXT 

1169 CURSET 226,14,3 

1174 FILL 164,2,18 

12GG PRINT "PASSAGE HUMERC: "5H 

1214 RETURH 


Fonction 
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17. TRI DE FICHIERS 





Il s'agit d'un processus important dans la manipulation de fichiers. 

Soit un ensemble de N nombres dans un ordre quelconque. Supposons- 
les placés dans une matrice A (N). Le but de l'algorithme est de ranger 
dans l'ordre croissant (ou décroissant) cet ensemble de nombres. 
Différentes méthodes existent; certaines sont beaucoup plus rapides que 
d’autres; le critère de sélection est naturellement le nombre N. Si ce 
nombre est grand, il faut utiliser une méthode sophistiquée et rapide; s’il 
est petit, une méthode plus simple convient. 

Le tri de mots est effectué de la même manière. En effet, pour un 
ordinateur, les lettres ont un ordre, identique à l’ordre alphabétique, 
défini par leur code ASCII. Ainsi : 


ARBALETE < ARBRE <FORET <ZOO 


Cependant, des problèmes apparaissent si minuscules et majuscules 
sont mélangées. On a : 


Majuscule < minuscule 


Ainsi : 


ue W<X <Y <Z<a<b< cc... 


Cela étant dit, les algorithmes restent les mêmes, dans la mesure où A (1) 
est remplacé par A$(1). 

L'opération la plus souvent effectuée dans ce genre de problèmes est 
l'échange de deux données. Certains micro-ordinateurs disposent pour 
cela de l'instruction SW AP. 
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On peut ranger les différentes méthodes de tri dans trois catégories 
principales : 


a. Tris par échange (ou par inversion) 


Le principe est de comparer et d’inverser (éventuellement) deux 
éléments de la liste suivant une règle fixe. 
Appartiennent à cette catégorie : 


— Ripple 
— Bubble 
— Shaker 
— Quick 


b. Tris par insertion 
Chaque élément est inséré à la bonne place. Citons : 


— Insertion directe 
— Shell 
— Shell-Metzner 


c. Tris par sélection 


On sélectionne et on classe l’élément le plus grand successivement dans 
plusieurs sous-listes, en arbres. Ce sont : 


— Sélection 


— Heap sort 


Nous verrons en détail ces différentes méthodes, en donnant des 
exemples commentés. Puis nous verrons l'intérêt d’avoir des listes 
classées : l’accès à un élément particulier en est grandement facilité. Nous 
aborderons enfin l’accès plus général en donnant simplement les 
algorithmes dans leurs grandes lignes. Il est en effet ici impossible de 
donner des programmes sans faire appel aux propriétés particulières du 
DOS de la machine et des types (séquentiels ou à accès direct) de fichiers 
qu'il peut traiter. 
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Méthodes de tri 


A. Tris par échange 


Ripple sort 


Il s’agit de faire (N—1) passes sur l’ensemble de la matrice, en 
effectuant la comparaison systématique de deux éléments consécutifs, en 
inversant la place de ceux-ci s’ils ne sont pas dans l’ordre. Le nombre de 
comparaisons effectuées est de l’ordre de N°72 (O(N72)). 

En particulier, si À (N) est déjà rangée, l’algorithme ne le «sait» pas et 
effectue donc le même nombre de comparaisons (ce qui est un 
inconvénient). Cependant, l’utilisation d’un test (lignes 130, 170, 200, 
210) permet d’éviter ce phénomène puisqu'il est ainsi possible de savoir 
si, au cours d’une passe, il y a eu permutation à un moment quelconque. 
Dans le cas contraire, c’est que le tri est terminé. 


Programme 





Exemple de fonctionnement 


Soit à classer une matrice de 15 éléments, dans l’ordre croissant 
d’abord, dans l’ordre décroissant ensuite. Le résultat obtenu est le 
suivant : 





En remplaçant (ligne 150) 
réalise un tri inverse : 





Attention 


Les exemples sont donnés à titre indicatif, et permettent d'interpréter 
la manière dont le classement évolue en fonction du temps. Ils ne sont pas 
significatifs de la vitesse à laquelle évolue le tri. 


Commentaires 


Étudions la manière dont fonctionne le programme. Pour N donné, on 
effectue N —1 comparaisons sur l’ensemble des éléments à trier. Au 
premier passage, chaque élément est comparé au suivant. Il y a échange si 
l’ordre n’est pas le bon. L’algorithme enrichit donc la tête de liste d’un 
seul élément (le plus grand parmi les éléments non triés) à chaque 
passage. On conçoit donc que le nombre (maximal) de comparaisons à 
effectuer est Nx(N — 1). 
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Cependant, étant sûr que la partie triée augmente d’un élément à 
chaque passage, on peut (ligne 210) réduire simultanément d’un élément 
la nouvelle liste à trier. 

Nous ajoutons un petit programme qui trie une liste contenant des 
noms et des chiffres, une fois par ordre alphabétique, et ensuite par ordre 
numérique. 


Exemple 

1 2 
DURAND 33 ABEL 35 31 ISOBAR 
DUPONT 32 DUPONT 32 32 DUPONT 
ABEL 35 —+ DURAND 33 33 DURAND 
ISOBAR 31 ISOBAR 31 35 ABEL 
LION 41 LION 41 41 LION 
Programme 


à REM RIPPLE PLUSIEURS CRITERES 

14 INPUT CHÔMERE He "0H 

CA CIN AGCHS.ACMS 

30 FOR I=z1 TO H 

44 PRIT MHONM AMI Ti & Mi: THPUT ACT 


PRINT "YALEUR ACT: æ Mi: THPUT A6 


12) 


HET 

IHPUT TRI 4%, À, FIH CASBSE 5 "OR 
IF Gé "AU AND Gé "EE" THEH EHG 
REM TRI RIFPLE 

F=H 

PisF: FOR Is1 T0 P1-1 

Life T4 1 ss LEA T 41 5 

IF Che" AND UgreAge TS THEM 164 
IF G$="E" AND Lr=ACT THEM 164 
ABC I + 1 SAT 1e ABC T LE 

A6 +1 SCT: ACT el PeT 

HEAT 

IF FisP THEH #64 

GOTO 138 


Un RL ii eo PR IR 


RD D RAM 


He He Hi Ha bé be De pe be He Lt OGC 


D JTE 
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284 FOR I=1 TO H 
CAS IF GS="E" THEM PRINT ACT. AI: GO 
TS 


LL 

214 FRIHTASET ACT 
2eG HET 

cn COTO 84 


Le paramètre P1 indique s’il y a eu inversion, sinon la liste est en 
ordre, le tri est terminé. Le paramètre P marque la position à partir de 


laquelle les éléments sont classés. 


Quand on trie une liste de chaînes A$(1) à l’aide d’un programme en 
BASIC il faut prendre en considération — surtout si le nombre N est très 
grand — le fait qu’il y a de temps en temps un ramassage de déchets 
(Garbage collection) qui prend beaucoup de temps. Il vaut donc mieux 
introduire un vecteur clé B(I) et ne pas inverser les variables chaînes. 


Voici le petit programme qui trie les indices au lieu des chaînes : 


à FREM RIFPFLE PLUSIEURS CRITEFES-ACCES F 
AR CLE 

14 IHFUT "HOMERE He "5H 

CA CIM A&CHS HCMI BCHS 

34 FÜR Isi TO H 

48 PRIHT "HOM ABOMITI MS = Mi: INPUT AI 


cn 


4 PRINT "YALELIR HE"iTi"u æ M; INPUT fé 
Jill] 
ÊA MESXT 
EG IMHPFUT TRI A$, 4. FIH CH°EB-F 5 "cé 
SA IF Gé< et AND Géé2"E" THEH EHC 
144 KEM TÉI KIFFLE 
124 FzH 
134 Pi=F: FOR Izi TO F-1 
LES Ve T4 1 5: be Age M: Le T41 
14G IF Gé="A" AND UfieA$é EC] THEN 16 
145 IF GS="E" AND UreAtTr THEH 164 
154 BC T+1 16 TI 5: BC I = 
1e AE I+A1 = TU: AE T ell:PaeT 
164 HEXT 
LFP IF Fi=F THEH 244 
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124 GOTO 134 
cAG4 FOR I=1 TO H 
PAS IF Géz"E" THEH FRIHT ACT HECECT 
:GOTOZZA 
214 PRINTAGEECT ACT 
BEA MEET 
£a GOT 54 
Remarque 


Un tri Ripple est très efficace sur une liste déjà ordonnée au départ. 


Tri Shaker 


On peut essayer d'améliorer le procédé du Ripple en commençant les 
comparaisons une fois par la gauche et immédiatement après par la 
droite. Le succès dudit Shaker sort n’est pourtant pas évident. 

Voici le programme : 

Les paramètres D et F limitent la sous-liste qui reste à trier. À noter 
une petite acrobatie de BASIC (Microsoft) autorisée par la méthode de 
stockage des variables D et F (valeurs initiales et finales des boucles) 
d'une boucle FOR... NEXT. 160. F=I ne modifie pas la valeur finale 
de la boucle en train. Vérifiez-la sur votre machine. 


14 FEM SHAKER RIFPLE 

INPUT "HOMERE H= "M 
CIM ACH 2 
FÜR I=1 TO H 
ACT: = IHTÉRHDé 41 1414845 
HET 

REM TRI SKAKER 

C=i:F=h 

TEST=6 

FOR Is TO F-1 


JOUR 
RER ENTER ENTER ENT 


mire pe Ii Le it Tri 
QD in 
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154 IFACI ec I+19 THEM 174 

16A Le ACT: ACT ACT 41 1: ACT 41 EL: TEST=I 
FT 

174 MEXT 

184 IF TEST =4 THEH 364 

183 TESTS 

194 FOR I=F TO D+1i STEF -1 

21 Ro Se 12 THEH Z54 

224 Us T ACT SAC I-1L HE ACI-1 IL TESTEL : 
(RE 

224 HEXT 


Pdf IF TEST=4 THEH 596 
254 GOTO1ZA 


249 FOR I=i TO H 
214 PRINT AI; 
324 MEXT 

334 STOP 


Tri Shaker : Exemple 





Commentaire 


Cette méthode de tri porte un nom qui indique bien la manière dont 
elle opère. Le premier passage est classique, du type Ripple. Le plus grand 
élément est amené en début de liste. Le second passage est un Ripple 
inverse. Partant du début de la liste, on ramène à la fin le plus petit. 
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Entre-temps, la liste en question a été réduite à la partie non triée, et cela 
se produira à chaque passage : 


partie triée 








Premier passage —+ [ 
Ripple direct t 


1 


L 











DS se partie triée 
Ripple imese — PSN 222 
. o 
Troisième passage —+ CD WW 





nouvelle valeur de F 


Ripple direct 


et ainsi de suite... 
On pourrait dire que F et D jouent le rôle de mémoires, conservant 
pour la prochaine étape du tri l’indice du tableau pour lequel le dernier 


échange a eu lieu. 
Bubble sort (tri à bulle) 


Le principe est le suivant : 


e Dans la liste des N nombres, on recherche le plus petit. Celui-ci est 
alors placé en tête (premier élément). 


e Dans la liste des (N—1) nombres restants, à partir du second, on 
recherche le plus petit. Il est placé en tête de cette nouvelle liste et 
devient en fait le deuxième élément de la liste triée. 


e Dans la liste des (N — 2) nombres restants, … etc. 


On conçoit que les éléments remontent, dans l’ordre du tri, comme des 
bulles dans un tube à essai, d’où le nom de la méthode. 


Programme 
Q FREM ne PR LE 
14 FEM CONSTRUCTION CE LA MATRICE 
11 REM -- DE ein où de ce ie 
£4 INPUT "HOMBRE H= "5H 
39 DIM ACH3 
43 FOR Izi TO MH 


— 258 - 


94 AC I =INTERHDC 1 3410008 

64 HEXT 

A0 RENE esse 

14 REM EUEBLE SORT 

1A1-HEM FS rennes 

124 FOR I=1 TO H-1 

159 FOR Jsl+1i TO H 

14 IF ACID=<AC I) THEH 164 
159 U=AC I 2: ACT SEA 4: A6 J'all 
169 MEXT J 

179 HEAT I 

184 STOF 

208 FOR I=1 TO H:PRINTACI 2; :HEXT 


Exemple de fonctionnement 


En reprenant la matrice précédente, on réalise ce programme. On peut 
constater que les éléments remontent dans l’ordre vers le haut : à la 
1ème passe, les I premiers éléments sont déjà rangés. 





Sn CE 
JET Un Re 


D CR nd 


mhmmbemebkmil 





BLELB ER LE 


Commentaire 


Le premier élément est comparé à tous les autres (et, le cas échéant, 
échangé avec ceux-ci). Le plus petit a donc obligatoirement pris sa place 
à la fin du premier passage. Le second élément est ensuite comparé à tous 
les autres et, par le même procédé, remplacé par le plus petit de la liste 
restante. Ce phénomène se poursuit jusqu’au classement final. 
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Premier passage 1 N 
ne se 
PR 





J varie de 2 à N le plus 
petit éléments examinés 
Deuxième passage 12 N 
I=2 
: N ces deux ————  —" 
J varie de 3 à N éléments sont ” 
maintenant classés éléments examinés 
Troisième passage 1 2 3 N 





ee SSI 54 Lu 
J varie de 4 à N ces éléments ———]“ —— 


sont classés éléments examinés 





Remarque 


On pourrait se demander quelle différence existe entre un tri Ripple et 
un tri Bubble. Celle-ci est subtile. Bubble compare successivement et 
systématiquement chaque élément de la liste à tous les autres. Ripple fait 
de même, mais sur deux éléments consécutifs. Cela fait qu’il se constitue 
à chaque passage des micro-listes triées, dans la mesure où un élément 
temporairement grand (devant ceux auxquels il a été comparé) remonte 
la liste jusqu’au moment où il en rencontre un plus grand que lui. Pour 
poursuivre l’analogie hydraulique, il y a remontée des petites «bulles »; 
celles-ci s'arrêtent dès qu’elles en ont rencontré une plus grosse. 

La méthode Bubble se justifie lorsque le tri n’a pas besoin d’être 
complet; on se contente, par exemple, de ne connaître que les trois 
«premiers » d’une liste. 


Quick sort (tri rapide) 


Il s’agit d’une méthode sophistiquée, très adaptée au traitement des 
fichiers importants. Dans ce cas, le temps de tri peut être nettement 
inférieur à celui que nécessitent les autres méthodes. 
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Programme 


19 


DONS ISPRIR IS 


GOTQ 19 

REM KXLEREEREEE RE ERELEEREELELELEEAAX 
REM # IHITIALISATION TABLE NOMBRES # 
REM ÆRAEENEEERERERERERERELEELETENEALE 
REM 

REA Ut AC EH SI: AS ISX 

RETURH 


FEM 
INPUT "NOMBRE M= "NM 


11 HMI=INTCS#LMOM 35 
12 DIM ACH° 
13 INPUT "LOHGEUR MICRO-LISTES 4=élF=<1 
INTER 

15 INPUT "TRI FINAL RIPFLE - INSERTION 
CRI) "5 TFS 


a 


FOR I=1 TO H 


ACID = INTÉRHMDC 10#10090 ) 
HEXT 

PRIHT "DEEUT TRI" 

KEM 


REM HHEAEAERERERELELEELELELELERE 
REM # IHITIALISATION VARIABLES # 
REM HERALELE RER EEE ERELELAEELEX 
REM 

DIM PLEC1>;M13 

PP=4 : BAS=1 : HAUT=H 

PYOT=G:MIL=A 


EN 
RE  ÆÉAEMERETELEREEARE 


198 REM *% CHOIX CU PIYOT # 

101 REM XXE EELEENEZREENEX 

192 FREM 

114 MIL2THTÉ EAS+HAUT 342 

124 IF ACBAS MACMIL) THEN U=RBAS : S=MIL:G 
JSUB? : J=A 

124 IF ACHAUT AC BAS: THEM U=BAS : S=HAUT 
:GOSUEF : I=1 

144 IF Is1 AH ACBASIEACMIL]Y THEN UEMIL 
:SeBAS: COSUE 

154 REM 
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191 REM HAE EEEEEEEEERELERERERENAEAE 
152 REM * CLASSEMENT CES ELEMENTS # 
52 FEM # AUTOUR CU PIVOT # 
194 REM ÆÉFELAEEEEREFEREEEEEELEEEEE 
1255 FREM 

164 PYOTS=AC ES 

174 PE=<EAS : FH=HAUT +1 

184 REPERT 

194 PE=FR+1 

244 UNTIL ACPFE A 3=PLOT 

214 FEPEAT 

224 FH=PH-1 

254 UMNTIL ACPHiéePYQT 

244 IF FES=FH THEH 274 


F2 Pit 
13 
er] 


Fit PQ Pi Fo Ti 
59 05 Q Qt Ti 
rs 1Q 1 


2 
re 


RD in 


= 
se 


D PO re QD O0 4 A M RÉ Po pa 10 LEO 


BR BR a Be 0 09 0 CO QD Ci LA 0 Lo Lo POP ND 


BE à 
Ln R 


8 R 
JD: 
RS 9 D DIR CN DR DR D DE: 


4sf 


LISPE : SeFH:GOSUET 

GOT 154 

LISBAS : SePH: GOSUE7 

FE M 

FEM HÉÉÉERERERERERERELEEEEEEAEX 
FENM # ORGANISATION CE LA FILE # 
FENM HER ARE HAE E A ACE AE LE LE ERENX 
REM 

SA re THEH 44 
IF CFH-EASéELF THEN 564 
FFEFF +1 

FLEC EG, FF 1=PH4t 

PLES 1. PF =HALIT 

HAUT =FH-1 

GOT 1 

IF CHAUT-FHS<ELF THEM SA 
EAS=PH+1 

GOTO11 

FEM 

IF CHAUT-FHi<ELF THEH 464 
FFEPP+1 

FLECG, PP <EAS 

PLES 1, PFi=PH-1 

EAS=FH+T 

GOTI 114 

IF CFH-EAS IÉELF THEH SG 


4 HAUT=FH-1 


GATO 114 
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RE 
RENM HÉHEAHAEEHE 


1e REM # DEPFILER *# 
3 REM KEFELELEEAEX 


RER ER EN D 
0 CO NN 9 00 D R 9 Dm: 


DS 9 DD fr Qi NT 


D'AIDE 
9 Firm: 


M 
Un R 
LAN EUX 


M 
LA | 


ID EN De DER: 


17 


to 15 


E", 


D A CP me 9 D ED NO Mt 9 CO D IQ C0 


CO NT TS ND DM DM MD DAMMPDE 


NS 0 19: 


© 


it 


ai 


FEN 

IF FP=4 THEH SEA 
EASS<PLES 4, PP 
HAUT=PLES1.FPû 

PF=PP-1 

GOTO 114 

IF TF&2="I" THEN Fam 

FRET 

RE HAE EEE EE RAR NA EM AUECE 
FEM %# TRI FIHAL KRIFFLE # 
RENM ÆHÉRERELELE RAA AA AAA EEE 
FE 

I:=M 

FRIHT'OEELUT RIFFLE" 

FEPFEAT 

TEST=A 

FOR .I=1 TO I-1 
IFACJN=< AE .I41 1 THEH €E@ 
ASE I: AC Vie UT re A 141 té 
TEST=1 

HET 

I=I-1 

UHTIL TEST=HA 

GOTO 544 

REM 

FEM HÉLEREE ER ERA AE HAN AAHAEE 
FENM %* TRI FIHAL INSERTION # 


REM 

PRIHT'CGERUT IHSERTIOH" 
FOR Je TA H 

IF Hé tesfAl.I-15 THEH FF 
LIsHé li: Get-1 

FEFEAT 

AG G+1 AG 51: AG el: GeG-1 
UHTIL Ge OR List Gi 
HET 


FOR Is TO M:PEIHT HI 3 :MNEXT 
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Commentaires 


Soit une liste de 256 éléments à trier, rangés dans A (256). 


Le programme débute ainsi : 


— Définition d’une table pile PLE (1,16) 
(N=256, NI=INT(3*xLN (N))= 16). 


— Définition d’un pointeur de pile PP (=0). 


— BAS et HAUT : pointeurs d’extrémité de partition à trier. 
(Au départ : BAS=0 et HAUT =N=256.) 


— PVOT : détermine la valeur du pivot (= 0). 
— MIL : pointe le milieu de la table à trier (= 128). 
— N : nombre d’éléments à classer (— 256). 


— PB, PH : pointeurs bas et haut. 


Le tri commence ligne 110 : 


| | 


BAS MIL HAUT 





— A(BAS) et A (MIL) sont comparés. Le plus grand des deux est 
placé en BAS; I=0. 


— A(HAUT) et A (BAS) sont comparés, le plus grand des deux est 
placé en HAUT; I=1. 


— A(BAS) et A (MIL) sont comparés si I= 1; en effet, dans ce cas, il 
n’est pas certain que A (BAS) > A (MIL). 


Prenons un exemple. 


Au départ : 9 14 6 I=0 
Après exécution de la ligne 120: 14 9 6 = 0 
Après exécution de la ligne 130 : 6 9 14 I=1 
Après exécution de la ligne 140 : 9 6 14 
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Le principe de la suite des opérations est de trouver un pivot, élément 
de la table ayant une valeur moyenne telle qu’il y ait par rapport à lui 
autant d’éléments plus petits que d’éléments plus grands. Dans notre 
exemple, le pivot est choisi parmi le premier, le dernier et le médian de la 
liste. Il est placé en tête de liste. 

On a donc : 


PVOT = A (BAS) 


e Tous les éléments sont classés par référence au pivot. Partant du bas 
de la liste, on remonte celle-ci jusqu’à trouver un élément plus grand 
que le pivot, d’indice PB. De même, partant du haut de la liste, on 
descend celle-ci jusqu’à trouver un élément plus petit que le pivot : 


Plus petits Plus grands 
que PVOT que PVOT 


pivot 
PB PH 
échange éven- 
tuel entre 
éléments 


Si PB>PH, toute la liste est répartie. Le programme peut se 
poursuivre. Sinon, les deux éléments A (PB) [plus grand que PIVOT] et 
A (PH) [plus petit que PIVOT] sont interchangés. Et le même cycle 
recommence. Cela continue jusqu’au moment où PB2>PH, c'est-à-dire 
quand tous les éléments plus petits que PVOT sont d’un côté de la liste et 
que tous les éléments plus grands que PVOT sont de l’autre côté. A ce 
moment (ligne 270), le PVOT est placé en séparateur entre les deux 
demi-listes. 


Plus petits que PVOT Plus grands que PVOT 
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Après chaque classement autour du pivot, quatre tests sont effectués. 
Au premier passage, la liste concernée est la liste primitive, scindée en 
deux autour du pivot. 


290 x  L,>L, —> 400 
L, <LP —+ 360 
Donc L,<L;, —+ mémoriser la position de L, et retrier L.. 


360 # L,<LP —+ 500 


Sinon retrier L;. 


400 # L,<LP —+ 460 


Sinon mémoriser L, et retrier L.. 


460 * L, <LP ——+ 500 


Sinon retrier L.. 


500 : Les divers tronçonnages de la liste primitive ont conduit à une 
liste de longueur inférieure à LP (4 LP<10). On vient maintenant 
récupérer les listes mises de côté à chaque étape. 

e Si la liste L, est plus longue que la liste L,, on renvoie en 400. 
e Si la liste L, contient moins de LP éléments, on renvoie en 360. 


e Le cas traité est donc : 
longueur de L, longueur de L, 
c'est-à-dire que la demi-liste la plus courte est traitée. On a alors : 


PP=PP+1=I1 
PLE(0,PP)=PLE(0,1)=PH+1 
(pointe le début de L,) 
PLE(1,PP)=PLE(1,1)=HAUT 
(pointe la fin de L,;) 
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Le programme renvoie alors en 110, pour effectuer un tri sur L;, 
comme il avait été effectué au départ sur L, +L, (choix d’un pivot dans 
L,, répartitions des éléments de L, autour de ce pivot). 

Cela se poursuit jusqu’au moment où la dernière répartition conduit à 
des micro-listes contenant moins de LP éléments. 


e En ligne 400, pour laquelle on a déjà L, >L,, la longueur de L, est 
testée : si elle est inférieure à LP, il y a renvoi en 460. Sinon : 


PP=PP+1=1 
PLE(0,PP)=PLE(0,1)=BAS 
PLE(1,PP)=PLE(1,1)=PH—1 


Il y a mémorisation des coordonnées de L, et retour en 110 pour 
traiter L, (choix d’un pivot dans L,, répartition des éléments de L, autour 
de ce pivot). 


e En ligne 460, si les micro-listes de L, sont d’une longueur inférieure à 
LP, il y a sortie de ce type de tri. Sinon, on recommence en 110. 


e En 500 : test de PP. Si PP = 0, c’est que toutes les micro-listes ont une 
longueur inférieure à LP; le tri est presque terminé. Sinon, il faut 
chercher quelle est la liste de longueur supérieure à LP et effectuer un 


nouveau traitement. 


e En ligne 580 nous disposons maintenant d’un ensemble de listes de LP 
éléments prétriées et prérangées. Un tri final Ripple ou Insertion 
effectue le dernier classement. Il est très efficace à ce moment-là. 


Nous préférons pour les listes importantes selon la proposition de 
D. Knuth LP=9 et un tri final Insertion. Le choix est légitimé par une 
étude mathématique. 


Prenons un exemple : 
Soit la liste de vingt nombres donnée ci-après. Le pivot est A (1) = 13. 
(C’est une coïncidence.) Il est donc déjà placé en tête de liste et l’on n’a 
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pas à effectuer ce transfert. Suivons le déroulement du programme pas à 
pas : 


1. x 84(2>13) est échangé avec 11 (<< 13) 


* 21 est échangé avec 2. 


Le passage est terminé; 2 est maintenant échangé avec le pivot 13. 
Nous avons maintenant affaire à deux listes; l’une L1, à gauche du 
pivot, de longueur inférieure à 4 est estimée prétriée. 

Elle ne sera plus touchée. 


2. On s'occupe alors de L2, s'étendant de A(5) à A(20). Comme 
précédemment, le pivot est choisi dans cette liste : 


20+5 
A) = au =82, A(5)=85, A(20)= 89 
d’où le pivot 85, placé en tête de L2. 
Alors : 


3. + 88 est échangé avec 84 


* Le pivot 85 est échangé avec la position pointée par PH, soit avec 
17: 


4. Nous avons maintenant deux nouvelles listes, que nous noterons L1 et 
L2, situées de part et d’autre du pivot. 
La plus longue, L1, est mémorisée, c’est-à-dire que ses indices de 
début et de fin sont placés dans une pile (tableau), qui permettra au 
programme par la suite de retrouver les tronçons de liste non triés 
pour pouvoir les traiter. 


5. Le pivot est maintenant choisi dans L2, entre 89, 88 et 99. Il vaut 
donc 89 et il est placé en tête de liste. 88 et 94 sont alors interchangés, 
puis 89 avec 88. Nous aboutissons à deux nouvelles listes, de part et 
d’autre du pivot, mais de longueurs inférieures à 4. Elles sont estimées 
prétriées. 
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6. La liste A (5)-A (14) est alors dépilée, le pivot est 24 (choix entre 24, 
14, 84). 


7. 42 est échangé avec 21, puis 17 est échangé avec 24. Nous 
aboutissons à deux listes, dont l’une (L1 de A(5) à A(9)) est de 
longueur inférieure à 4. Elle est abandonnée. 

L2 (de A (10) à A (14)) est triée : pivot 82. Un seul échange a alors 
lieu : 82< — >73. 


8. Le tri Quick est fini. Un tri Ripple conclut. 




















(o) 

13 74 4 14 97 © 70 27 73 17 11 4 99 24 55 
13 74 42 14 97 2 70 37 73 17 11 4 97 24 £&S 
HAUT= 70 FE= © FH= 0 

O) FPIVOT= 13 

13 4 84 21 14 97 79 73 94 97 24 

13 4 11 21 14 77 79 73 74 9% 24 

13 4 11 2 14 77 70 73 4 97 24 

2 4 11 13 14 97 70 73 94 97 24 

IN LE TRI 

©) 
LI est triee: 

2 4 11 42 14 97 21 70 SZ 73 35 17 4 94 979 24 87 
HAUIT= :0Q FE= © FH= 4 

2 4 11 13 70 73 34 33 9 

ë 4 ll, 13 70 73 94 97 97 

& 4 11 13 70 73 94 97 77 
FIN LE TRI 

@ 
Li est memorisee 
LI est triees: 14 

£ 4 11 13 14 24 21 70 SZ 73 34 2S 7 94 99 
HALIT= 20 FE= 14  FH= 15 

(& 

(5) FIVOT= 27 

= 4 11 13 17 12 427 14 21 70 73 97 99 
É 4 11 13 17 15 47 14 Z4 Z1 70 73 97 77 
FIN DE TRI 

Fr 

@) 
La liste = 5 —- 14 est depilee. 

> 4 ii 13 T4 15 47 14 17: 21 79 82 73 34 35 BE 57 74 97 77 
HALIT= 14 BAS= ©  PE= 15 FH= 17 

7 

® FIVOT= 4 

4 11 13 24 1$ 21 14 17 42 70 73 7 27 97 
4 11 013 17 1% 21 14 24 42 70 73 97 99 





FIN DE TRI 
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5 


LI ect triee: 10 - 14 

Z 4 11 17 17 15 21 14 24 &% 70 42 73 &4 
HALT= 14 BAZ= 10 FE= 10 FH= 9 

(8) FIVOT- &2 

Z 4 11 173 17 13 Z1 14 24 73 70 42 82 64 
FIN DE TRI 


TRI TERMINE 


B. Tris par insertion 


Insertion directe 


85 SE 857 94 97 97 
ESS SE #7 94 97 97 


Soit une liste de N éléments à trier. On inverse les éléments A (1) et 
A (2), si nécessaire, et on range après les éléments A (3), A (4), …., A(N) 
en commençant par A(3) et en comparant aux précédents. Voici le 


programme : 


REM TRI 


LHELIT 





li = 


FRTHT 










4 


LITHEN 1 





REFEAT 
FE G+1 ve 
LHTIL 

NE 
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Tri par insertion : Exemple 





M LA 






ARR E 







pi pe he fee pou 
ini 


pa 
he 
£ 


Tri Shell 


Cette méthode rappelle le tri à bulles. La liste est scindée en deux 
parties (à peu près égales). Les premiers éléments sont comparés aux 
derniers. Puis la liste est scindée en quatre. Même opération, et ainsi de 
suite. Le pas de travail est ainsi de plus en plus petit, ce qui conduit à une 
diminution dynamique du nombre de comparaisons. 


Programme 








Exemple 





15 4" rd ; 
15 4 F1: 
15 4 4 
15 + 4 
s 4 4 
15 4 da 
4: 
d' 
4° 





Commentaires 


ns idée tu ; Fou" 


ime Litier Le False: 





tr au début LE 


Première Role 


U=d 





Un chanienment 


En Le 








soit successivement : 


PSE M 
un Chansement 






re 


AÉRER LE Lt 





min un 


= 4 


me ire iii 


pes Li ai LR 


in 


{es ji espe pes 


4 


SE ME 


#7 
à aucun changement 





# Auatriéme bon 


Det ri terils 


+ 
. 





Tri Shell-Metzner 


I1 débute comme le tri Shell (scindage en deux moitiés de la liste), mais 
à la fin du premier passage, chaque demi-liste est triée séparément. 


Programme 


REH TÉTI SHELL METZHER 
IHPFUT" HOMÉÈRE Hz "iH 
DIM ACM3 

FOR I =1 TO H 

HT = THTÉ BND 1 166 
FRIHT ACT 5: 

MEXT 

FRIHT 

REM SHELL METEHER 

D= IHTCMHAZ 

1: E=H-f 

=. 


ni 


et pb Hs pe pt js IT y LA A D CD fitr 
D LA Re 0 re CT SG CT ST QD 


IS Qt 9 1Q 
ei ÎLE 

1 

+ 
© 


Ti H 
ZE 
{+ 
il 
T 
LS 
m1 
— 
L 
ol 
+ 
H 
iQ 
1 
mn 
— 
2 
Fr 
£t 
13 


CON ACC = AE ACE DS 


su 


nd 
DA 
EU 


F Cè=1 THEH 
Je.t+1 

IF JI=£E THEH 
C=THTED4Z 
IF Dés THEH 
FOR I=1 TO H 
PRINT ACT: 
HE 


Fi 
bb 
L 
ii 


Pre 


Fr rares 
Je 19 D 1 —J 
SD oo 19 


Fes 
[RS | 
Len 


Ji 


F3 Po Fe Pi 
NiNBR ET 
DNERNELT 


Nous supposons, comme dans les cas précédents, que la matrice à trier 
est A(15). 





C. Tris par sélection 


On cherche l'élément le plus grand dans la liste comprenant F éléments 
et on le range à la tête, F=N, N—1, …, 1. 


Commentaire 


Dans un premier temps, toute la liste est examinée. Le plus grand 
élément, U, est amené en début de liste. Au passage, l’élément qui était à 
cet emplacement est amené à l'endroit où a eu lieu la dernière 
comparaison positive. La liste est ensuite réduite d’un élément et le 
cheminement du programme reprend sur cette liste raccourcie. 


Programme 








Exemple 


Si 





Exemple détaillé. 
Reprenons la liste utilisée précédemment : 


15 6 21 4 16 17 5 10 

Au premier passage, 21 est échangé avec 10 : 
15 6 10 4 16 17 5 21 

Au deuxième passage, 17 est échangé avec 5 : 
15 6 10 4 16 5 17 21 
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Puis successivement : 


15 6 10 4 5 16 17 21 
5 6 10 4 15 16 17 21 


5 6 4 10 15 16 17 21 
5 4 6 10 15 16 17 21 
4 5 6 10 15 16 17 21 


La méthode est particulièrement bien expliquée dans le livre de 
D. Knuth : Sorting and Searching (pp. 145-149). 

Notre programme commenté se fonde sur l’algorithme H (pp. 146- 
147). L'idée consiste à créer des sous-arbres ordonnés... 


2 3 
: es ail D — 
Dr 7 4 FR: 
AD ue 
fils fils 
16 I=1,2,..,17 


ce qui veut dire : 
I | I 
A(INT ())>410 si 1<G< INT ()<1<) 


… et à garder cette structure favorable en remontant les grands éléments 
à la tête. 


Sous-arbre rangé 


(a12) 
(17) (318) A (4), A(8), A(9), A(16), A(17) 


Père 


fils fils aîné 
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Algorithme H 


Soit une liste de N éléments à trier, rangés dans A (1), …., A (N). A la 
fin on aura A(1)S A(2)<, …, SA(N). 
Le programme débute ainsi : 


e définition de deux pointeurs d’extrémité : 
G pointe sur la gauche de la liste, 
D pointe sur la droite, 
N 
au départ : G=INT (>) +1, D=N (ligne 100) 


e la variable A sert de mémoire. 


Le tri commence ligne 200 du programme. 

Dans un premier temps, on diminue G, et ensuite, arrivé à G=1, D 
jusqu'à D =1. 

A ce moment-là, le tri est terminé si on range l’élément A en bas de 
liste : A(1)= A. 


Pendant la première phase (G>1), on crée des sous-arbres prétriés, 
c'est-à-dire : 


A (nr ()) > A(I)(x) 


si : 
I 
G< INT ()<1<x 


A la fin de ce premier cycle, tous les sous-arbres sont ordonnés et (x) 
reste valable avec G= 1. L'élément le plus grand est évidemment en tête 
de liste (= A(1)). 


On peut donc réduire la liste en posant A=A(D), (D=N), 
A(D)=A(1), D=D—1 (ligne 210). Dans la liste restante, qui com- 
prend D éléments, c’est-à-dire un élément de moins, on fait remonter en 
tête le plus grand élément, qui sera ensuite placé sur A (D), (D=N—1), 
après avoir posé A = A (D), etc. 
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A la fin, on arrive à D=1 et l’élément le plus petit devient A (1). 
Comment remonte-t-on les éléments? 


En première phase, on a A = A (G) et ensuite A = A (D), dont il faut 
trouver la bonne place dans le sous-arbre approprié en préservant le 


préordre. 


Dans ce but, on pose d’abord J =G (ligne 300) et ensuite (ligne 400) : 


I=J et 
J=2x*]J si A(2*x]J) >A(2*xJ +1) ou 
J=2*xJ+1 si A(2*J)< A (2*J +1) (ligne 500) 


(A (J) devient ainsi le fils aîné du père A (1).) 


Si A> A (J) (ligne 600), on a trouvé la bonne place pour A (l’ordre est 
respecté) et le programme nous renvoie ligne 200 après avoir posé 


A(1)=A (ligne 800). 
Sinon, on descend par 1=J, J—2*J ou J= 2x] +1 etc. et on reprend 


la boucle. 


Prenons un exemple (N = 5): 
Soit la liste de cinq nombres donnée ci-après : 


97 10 1 77 S5S 


OÙ : 


ESC 


10 
TT SS5 


On a G=2, D=5, A=A(G)=10 et on échange 10 avec 77 : 


RE 


77 
107 55 


G=1, rien ne se passe, la première phase est terminée, l’arbre est 


ordonné. 
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On réduit la liste en posant : 
A=55, A(5)=97, D=4 
I=1, J=2 
A(2) —+ A1) 
et on reprend : 


1=2, J=—4 


On a A>10—=A(4), donc À —+> A (2). 
Il en résulte : 


77 
es  — 


5 
107 
On pose A=10, A(4)=77, D=3 et 10 reste en bas: 
55 
Maintenant : 
A=1, A(3)=55, D=2 
et finalement : 


A(2)=10 et A(1)=1 


19 REM HEAFSORT C'APRES D. KHUTH 

24 INPUT UHOMERE Hz ":H 

34 CIM ACH3 

44 FOR Ii TO H 

26 Hole INTCRHDE 1 1#1@G4 5 

A OHEZAT 

146 GEIHTIH-2 541: DzH 

119 REM CEEUT . G GAUCHE , D CROITE 
SGA IF Gi THEN GeG-1: AsACG OS: GTI 394 
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214 A=AC CD 1: AE Cie 1e Dei: IF Dei THE 
H 461 3=fA: GÜTOSAG 

224 REM DIMIMUER G ET APRES D JUS" A 
Les FIH Ci 


216 


J=G 
FENM ACIHTÉK EN SACKU, GÉINTEK SRE 
PREPARATION FOUR TRI SOLUSARERES 
Iz' JISF#I 

IF ID THEH 544 

IF IC THEH 644 

IF .J:D THEN #44 


A FEM TRI SOUCARERE 


IF AC ACL THEH JIel+1 
RENM L'AUTRE FILS CI PERE ACT FLUS 


GANG 


Es 


LD 0 Qi 0 CO J —J me 
RS EE UE NE 
EURRURENELR ET EL EUX) mi 


IF A%= AG. THEH 564 

FENM COMPARAISON PLACE TROUVER 
ACT I=AC I: GOTO di 

EN SOUSTET FAS FIMI, OH REMONTE 
AS Ie: GOTO 244 

RENM SOUSTRI TERMINE 

FOR I=1 TO H 

FEINT ACT: 

HEAT 
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DEA 15 ZE ZE 17 E94 SEé 444 44 LISTE 
.GDA 





Se 1 sf Séé A3 SE 





Ai 15 DÉS LE O1ET 94 ddé qua GEI 
454 15 SE 1Er 294 46 AB SE 
454 15 Re A ÜE AIT 1 
154 15 Se 17 £eSsg AUS EI 









nr y 





A 5 44 15 1 2e 15, 
1 1. 
4 à Hé id 15 St Se À 
1: 
4 ë LE 
Î 
4 LE 12% 
1 


AA du SUR du ré F1 
jo lat Sie Sd €7e 713 





n 34 SEE Sd 676 F1 
! fat Si DR LOT NE 





Ér J LISTE PRECEDENTE N 
G D A 


Il 
22 





ie pee 





1 En 

ETS Lie S$ 
+. 4 56 

Pr GS 1 18 
L 24, 6 

LE li ete in 
d CE 

dE TS Fo 6 
Î l î 

1 14 5 Fr 37 
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Exemples, tableau comparatif 





Les tris Quick, Shell-Metzner et Heap sort sont particulièrement 
efficaces et recommandables. 
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18. ACCÈS A UN ÉLÉMENT DE FICHIER 





Le traitement de fichiers s'effectue souvent à l’aide d’un lecteur de 
disquettes, associé à son DOS (Disk Operating System). Celui-ci permet 
en général de créer des fichiers à accès séquentiel (les rubriques sont 
disposées les unes à la suite des autres) et des fichiers à accès direct (à 
chaque rubrique est associé un numéro d’enregistrement). Nous enten- 
dons par rubrique un ensemble de données ayant un point commun. Pour 
un annuaire, par exemple, une rubrique contient : nom, prénom, adresse, 
numéro de téléphone. Si le fichier est important, la recherche d’une 
rubrique peut se faire : 


e séquentiellement, pour un fichier à accès séquentiel : il faut lire dans 
l’ordre toutes les fiches jusqu’au moment où l’on trouve celle que l’on 
cherche. Une telle méthode est très pénalisante en temps si le fichier 
est très important; 


e directement, pour un fichier à accès direct, si l’on connaît le numéro de 
la rubrique cherchée, ce qui n’est pas le cas en général, car à priori rien 
n’associe le nom « DUPONT » à la rubrique n° 15! 


Nous proposons deux méthodes permettant donc d’aboutir à un tel 
résultat. 


ACCÈS PAR CLÉ 


On constitue en mémoire centrale une table d’index. Celle-ci est une 
matrice de chaînes; l’indice de chaque variable est justement le numéro 
de rubrique du fichier à accès direct sur disquette. 
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A$(1)= «DUPONT » 
A$(2)= «DURAND » 
A$(3)=«DURANT » 


A$ (N) = « ZORGON » 


Pour un fichier figé, une telle matrice est définitive. La recherche du 
numéro de rubrique est alors aisée : 


LA INPUT "HO CHERCHES 7": 4% 
24 FOR Dei TO H 

Ba IF ASCLUSAS THEN Hell 

44 HEXT L 

Si FÉINT"LE HUMERO CHERCHE EST : “aim 





Il suffit maintenant d’aller chercher sur la disquette la rubrique M pour 
pouvoir l’exploiter. 

Le traitement d’un fichier dynamique est plus délicat; il faut pouvoir 
faire des suppressions et des rajouts. 


e Pour un ajout, il faut transférer le fichier en mémoire centrale, le 
compléter avec la nouvelle rubrique puis le trier, on peut alors créer la 
nouvelle table d’index. 


e Pour un retrait, l'opération est plus simple, puisqu'il suffit alors de 
décaler vers le haut l'ensemble des rubriques suivant celle qui a été 
supprimée. La table d’index est alors construite. 


HASH-CODING 


La difficulté, comme nous l’avons vu précédemment, est en général 
d'associer le numéro de rubrique à la clé, c’est-à-dire la chaîne de 
caractères correspondante. 


mi OA = 


On peut donc associer à chaque clé un numéro (qui sera celui de la 
rubrique concernée), déduit des propriétés de la clé elle-même. On 
pourrait ainsi faire correspondre à chaque mot un nombre qui serait, par 
exemple, la somme des positions de chacune de ses lettres dans 
l’alphabet : 


D—+ 4 D 4 
Des 2% À A 
RS "TS (D: 16 
A 1 0 4 
N —+> 14 N 14 
Dis 4, “1 30 

62 90 


Rubrique 90 : DUPONT Jacques, … 
Rubrique 62 : DURAND Paul, … 


Deux problèmes surgissent alors : 


e Il faut prévoir que le nombre N trouvé dans ces conditions puisse être 
compris entre 1 et le nombre maximal K de rubriques du fichier; on 
est d’ailleurs obligé de surdimensionner celui-ci, puisqu'il n’y a aucune 
chance que, à K mots, correspondent K nombres tous différents et 
compris entre 1 et K! 


e En outre, et c’est une conséquence, il n’est pas impossible (il est même 
très probable si le fichier est important) qu’à deux clés différentes 
corresponde le même code. On dit alors qu’il y a collision. Ce risque 
est d’autant plus réduit que le nombre de rubriques du fichier est plus 
grand que le nombre réel de rubriques. Cependant, une telle 
éventualité n’est jamais à écarter. 


La méthode de gestion des collisions peut alors être la suivante : on 
associe à chaque clé, dans le fichier, un numéro. Celui-ci sera égal à zéro 
s’il n’y a pas de collision, ou bien au numéro de rubrique de la clé portant 
le même code. 
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Exemple 


On recherche Dupont Jacques, sachant que le fichier contient en sus 
Dupont Paul et Dupont Pierre. Le code n’est calculé que sur le nom, soit 
par exemple 90. 


ce n’est pas Dupont Jacques 


211 228 Dupont Paul 


ce n'est pas Dupont Jacques 


228 Dupont Jacques 


Ainsi, pour un fichier de 200 rubriques, on peut affecter 250 places, les 
indices 201 à 250 correspondant alors au fichier des collisions. On se 
reportera à la référence citée précédemment pour plus de détails. 
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A 


Accès à un élément de fichier, 283 
Accès par clé, 283 

Algorithme de Cooley-Tukey, 167 
Algorithme de Routh, 111 


Approximation en moyenne, meilleure, 


137 
Archimède, la loi d’, 28 


B 


Boucle à verrouillage de phase, 109 
Bubble sort, 258 


C 


Chaînage de programmes, 65 
Charge d’un condensateur, 29 
Circuit RC, 19, 29 

Corrélation, coefficient de, 46, 47 


D 


Déphaseur RC, 103 

Dérivée numérique, 62, 63 
Déterminant, 94 

Développement asymptotique, 55 
Diagramme de Bode, 12 


E 


Équation 
linéaire, système d’, 94, 185 
non linéaire, 88, 179 
Équation différentielle 
de Liénard, 181 
linéaire du premier ordre, 19 
système homogène, 185 
linéaire d’ordre n, 117 
système du premier ordre, 195 


INDEX 


du second ordre, 179 

aux dérivées partielles, 198 

paraboliques-hyprboliques-elliptiques, 
198, 208, 226 


Erreur locale 184 


F 


Filtrage, 12, 26, 29 
Filtres passe bas, passe haut, passe bande, 


12 


Fonction 


de Bessel, première, seconde espèce, 52, 
55, 62 

erreur, 56, 57 

gamma, 55, 58 

de transfert, 12 

de Spline, 154 


Formule 


de Cramer, 102 
de Newton, 89 
de Runge-Kutta, 184 


Fourier, 


transformée de, 164 
série de, 164 


Fréquence de coupure, 104, 126 
Fraction continue, 57 


H 


Hash-coding, 284 
Haute résolution, 9 
Heap sort, 276 

Horner, schéma de, 121 


I 


Intégration, 22, 71 
Interpolation, 153, 164 


L 


Lignes équipotentielles, 33 
Lissage, 154 
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M 


Matrice 


complémentaire, inverse, de Hilbert, 
94, 105 

triagonale, 139 

creuse, 140 


Méthode 


de Bairstow, 117, 122 

de Crank-Nicholson, 203 

de Fadeev, 96 

de Fehlberg, 185 

de Gauss-Seidel, 149, 230 

de Liebmann, SOR, 230, 239 
de Newton, 121 

de Romberg, 73 

du trapèze, 22, 73 

de Van de Vel-King, 89 


Module, 13 
Mouvement ponctuel newtonien, 37 


N 


solution exacte, approchée, 38, 40 


Nombres complexes, 169 


O 


Ordre de la méthode, 63 
Orthonormalisation, 132, 133 


P 


Papillon, 168 

Passe bas, haut, bande, 12 
Pendule pesant, 71 

Phase, 13 

Polynüine de phase, 166 
Poisson, équation de, 226 
Polynôme caractéristique, 94 
Pont 


de Wien, 180 


de Wheastone, 130 
Procédé d’orthonormalisation 
Schmidt, 131 


Q 
Quick sort, 260 


R 


Racines (complexes), 117 
Règles de l'Hôpital, 62 
Régression linéaire, 46 
Résidu, 231 

Ripple sort, 252 


S 


Schéma de Horner, 119 

Série de Fourier, 164 
Shell-Metzner, tri de, 273 
Signal modulé en fréquence, 52 
SOR, méthode de, 230, 239 
Spectre, 164, 175 

Stabilité d’un polynôme, 109 


T 


Trace, 95 

Traitement de fichier, 283 
Transmittance, 12 
Transformée de Fourier, 164 
Tri, 250 


V 


Vibration d’une corde, 208 


Z 


Zéro d’un polynôme, 117 


— 288 - 


de 


Cet ouvrage propose une collection 
d'algorithmes numériques, destinés à être 
employés sur micro-ordinateur dans le 

but de résoudre nombre de problèmes 
scientifiques (mathématiques, physique, 
électronique). Les programmes de ce livre 
sont écrits en BASIC Microsoft et ne font 
appel à aucune propriété particulière d'une 
machine. 
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